Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2230)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2232)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adjustments to new topography concept
 ! 
 ! Former revisions:
@@ -139,10 +139,10 @@
     USE arrays_3d,                                                             &
         ONLY:  dzu, dzw, e, e_p, nr, pt, pt_p, q, q_p, qr, u, u_p, v, v_p,     &
-               w, w_p, zu, zw, z0
+               w, w_p, zu, zw
 #else
    USE arrays_3d,                                                              &
         ONLY:  dzu, dzw, e, e_p, e_1, e_2, nr, nr_2, nr_p, pt, pt_p, pt_1,     &
                pt_2, q, q_p, q_1, q_2, qr, qr_2, s, s_2, u, u_p, u_1, u_2, v,  &
-               v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw, z0
+               v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
 #endif
 
@@ -151,5 +151,5 @@
                message_string, microphysics_seifert, nest_bound_l, nest_bound_r,&
                nest_bound_s, nest_bound_n, nest_domain, neutral, passive_scalar,& 
-               simulated_time, topography, volume_flow
+               roughness_length, simulated_time, topography, volume_flow
 
     USE cpulog,                                                                 &
@@ -161,6 +161,5 @@
     USE indices,                                                                &
         ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg,  &
-               nysv, nz, nzb, nzb_s_inner, nzb_u_inner, nzb_u_outer,            &
-               nzb_v_inner, nzb_v_outer, nzb_w_inner, nzb_w_outer, nzt
+               nysv, nz, nzb, nzb_max, nzt, wall_flags_0
 
     USE kinds
@@ -201,4 +200,7 @@
 
 #endif
+
+    USE surface_mod,                                                            &
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -1310,4 +1312,5 @@
 
        INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i.
+       INTEGER(iwp) ::  end_index    !:  End index of present surface data type
        INTEGER(iwp) ::  i            !:
        INTEGER(iwp) ::  icorr        !:
@@ -1319,13 +1322,25 @@
        INTEGER(iwp) ::  jw           !:
        INTEGER(iwp) ::  k            !:
+       INTEGER(iwp) ::  k_wall_u_ji    !:
+       INTEGER(iwp) ::  k_wall_u_ji_p  !:
+       INTEGER(iwp) ::  k_wall_u_ji_m  !:
+       INTEGER(iwp) ::  k_wall_v_ji    !:
+       INTEGER(iwp) ::  k_wall_v_ji_p  !:
+       INTEGER(iwp) ::  k_wall_v_ji_m  !:
+       INTEGER(iwp) ::  k_wall_w_ji    !:
+       INTEGER(iwp) ::  k_wall_w_ji_p  !:
+       INTEGER(iwp) ::  k_wall_w_ji_m  !:
        INTEGER(iwp) ::  kb           !:
        INTEGER(iwp) ::  kcorr        !:
        INTEGER(iwp) ::  lc           !:
+       INTEGER(iwp) ::  m            !: Running index for surface data type
        INTEGER(iwp) ::  ni           !:
        INTEGER(iwp) ::  nj           !:
        INTEGER(iwp) ::  nk           !:
        INTEGER(iwp) ::  nzt_topo_max !:
+       INTEGER(iwp) ::  start_index  !:  Start index of present surface data type
        INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
 
+       REAL(wp)     ::  z0_topo      !:  roughness at vertical walls
        REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
 
@@ -1339,6 +1354,37 @@
           DO  i = nxl-1, nxl
              DO  j = nys, nyn
-                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),   &
-                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
+!
+!--             Concept need to be reconsidered for 3D-topography
+!--             Determine largest topography index on scalar grid
+                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on u grid
+                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on v grid
+                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on w grid
+                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
              ENDDO
           ENDDO
@@ -1350,6 +1396,37 @@
           i = nxr + 1
           DO  j = nys, nyn
-             nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),      &
-                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
+!
+!--             Concept need to be reconsidered for 3D-topography
+!--             Determine largest topography index on scalar grid
+                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on u grid
+                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on v grid
+                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on w grid
+                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
           ENDDO
           nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
@@ -1360,6 +1437,37 @@
           DO  j = nys-1, nys
              DO  i = nxl, nxr
-                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),   &
-                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
+!
+!--             Concept need to be reconsidered for 3D-topography
+!--             Determine largest topography index on scalar grid
+                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on u grid
+                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on v grid
+                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on w grid
+                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
              ENDDO
           ENDDO
@@ -1371,6 +1479,37 @@
           j = nyn + 1
           DO  i = nxl, nxr
-             nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),      &
-                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
+!
+!--             Concept need to be reconsidered for 3D-topography
+!--             Determine largest topography index on scalar grid
+                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on u grid
+                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on v grid
+                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
+!
+!--             Determine largest topography index on w grid
+                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
+                  MAXLOC(                                                      &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )  &
+                               ), DIM = 1                                      &
+                        ) - 1          )
           ENDDO
           nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
@@ -1422,9 +1561,50 @@
 !--          Left boundary for u
              i   = 0
-             kb  = nzb_u_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_u_l(1,k,j) = lc
              logc_ratio_u_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
@@ -1433,9 +1613,50 @@
 !--          Left boundary for v
              i   = -1
-             kb  =  nzb_v_inner(j,i)
-             k   =  kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_v_l(1,k,j) = lc
              logc_ratio_v_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
@@ -1469,9 +1690,51 @@
 !--          Right boundary for u
              i   = nxr + 1
-             kb  = nzb_u_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
+
              logc_u_r(1,k,j) = lc
              logc_ratio_u_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
@@ -1480,9 +1743,50 @@
 !--          Right boundary for v
              i   = nxr + 1
-             kb  = nzb_v_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_v_r(1,k,j) = lc
              logc_ratio_v_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
@@ -1516,9 +1820,50 @@
 !--          South boundary for u
              j   = -1
-             kb  =  nzb_u_inner(j,i)
-             k   =  kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_u_s(1,k,i) = lc
              logc_ratio_u_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
@@ -1527,9 +1872,50 @@
 !--          South boundary for v
              j   = 0
-             kb  = nzb_v_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_v_s(1,k,i) = lc
              logc_ratio_v_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
@@ -1563,9 +1949,50 @@
 !--          North boundary for u
              j   = nyn + 1
-             kb  = nzb_u_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_u_n(1,k,i) = lc
              logc_ratio_u_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
@@ -1574,9 +2001,50 @@
 !--          North boundary for v
              j   = nyn + 1
-             kb  = nzb_v_inner(j,i)
-             k   = kb + 1
-             wall_index = kb
-             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,      &
-                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--          For loglaw correction the roughness z0 is required. z0, however, 
+!--          is part of the surfacetypes now, so call subroutine according
+!--          to the present surface tpye. 
+!--          Default surface type
+             IF ( surf_def_h(0)%start_index(j,i) <=                            &
+                  surf_def_h(0)%end_index(j,i) )  THEN
+                start_index = surf_def_h(0)%start_index(j,i)
+                end_index   = surf_def_h(0)%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_def_h(0)%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_def_h(0)%z0(m), &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Natural surface type
+             ELSEIF ( surf_lsm_h%start_index(j,i) <=                           &
+                      surf_lsm_h%end_index(j,i) )  THEN
+                start_index = surf_lsm_h%start_index(j,i)
+                end_index   = surf_lsm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_lsm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_lsm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+!
+!--          Urban surface type
+             ELSEIF ( surf_usm_h%start_index(j,i) <=                           &
+                      surf_usm_h%end_index(j,i) )  THEN
+                start_index = surf_usm_h%start_index(j,i)
+                end_index = surf_usm_h%end_index(j,i)
+                DO  m = start_index, end_index
+                   k          = surf_usm_h%k(m)
+                   wall_index = k - 1 
+                   kb         = k - 1
+                   CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,  &
+                            j, inc, wall_index, surf_usm_h%z0(m),    &
+                            kb, direction, ncorr )
+                ENDDO
+             ENDIF
              logc_v_n(1,k,i) = lc
              logc_ratio_v_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
@@ -1590,8 +2058,20 @@
 !--    Then vertical walls and corners if necessary
        IF ( topography /= 'flat' )  THEN
+!
+!--       Workaround, set z0 at vertical surfaces simply to the given roughness
+!--       lenth, which is required to determine the logarithmic correction 
+!--       factors at the child boundaries, which are at the ghost-points. 
+!--       The surface data type for vertical surfaces, however, is not defined 
+!--       at ghost-points, so that no z0 can be retrieved at this point. 
+!--       Maybe, revise this later and define vertical surface datattype also
+!--       at ghost-points.
+          z0_topo = roughness_length
 
           kb = 0       ! kb is not used when direction > 1
 !        
 !--       Left boundary
+
+!
+!--       Are loglaw-correction parameters also calculated inside topo? 
           IF ( nest_bound_l )  THEN
 
@@ -1599,15 +2079,48 @@
 
              DO  j = nys, nyn
+                k_wall_u_ji   = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,0), 26 ) &
+                                 ), DIM = 1                                    &
+                                      ) - 1
+                k_wall_u_ji_p = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j+1,0), 26 )&
+                                 ), DIM = 1                                    &
+                                      ) - 1
+                k_wall_u_ji_m = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j-1,0), 26 )&
+                                 ), DIM = 1                                    &
+                                      ) - 1
+
+                k_wall_w_ji   = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,-1), 28 )&
+                                 ), DIM = 1                                    &
+                                      ) - 1
+                k_wall_w_ji_p = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j+1,-1), 28 )&
+                                 ), DIM = 1                                    &
+                                      ) - 1
+                k_wall_w_ji_m = MAXLOC(                                        &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j-1,-1), 28 )&
+                                 ), DIM = 1                                    &
+                                      ) - 1
+
                 DO  k = nzb, nzt_topo_nestbc_l
+
+                   i            = 0
 !
 !--                Wall for u on the south side, but not on the north side
-                   i  = 0
-                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.         &
-                        ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )            &
+                   IF ( ( k_wall_u_ji > k_wall_u_ji_p ) .AND.                  &
+                        ( k_wall_u_ji == k_wall_u_ji_m ) )                     &
                    THEN
                       inc        =  1
                       wall_index =  j
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1620,11 +2133,10 @@
 !
 !--                Wall for u on the north side, but not on the south side
-                   i  = 0
-                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND.         &
-                        ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
+                   IF ( ( k_wall_u_ji > k_wall_u_ji_m ) .AND.                  &
+                        ( k_wall_u_ji == k_wall_u_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  j + 1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1635,13 +2147,14 @@
                    ENDIF
 
+                   i  = -1
 !
 !--                Wall for w on the south side, but not on the north side.
-                   i  = -1
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) )  THEN
+
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_p )  .AND.                 & 
+                        ( k_wall_w_ji == k_wall_w_ji_m ) )  THEN    
                       inc        =  1
                       wall_index =  j
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1654,11 +2167,10 @@
 !
 !--                Wall for w on the north side, but not on the south side.
-                   i  = -1
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_m )  .AND.                 &
+                        ( k_wall_w_ji == k_wall_w_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  j+1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1682,13 +2194,45 @@
 
              DO  j = nys, nyn
+
+                k_wall_u_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j,i), 26 ) &
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_u_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j+1,i), 26 )&
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_u_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j-1,i), 26 )&
+                                     ), DIM = 1                                &
+                                        ) - 1
+
+                k_wall_w_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j,i), 28 ) &
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_w_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j+1,i), 28 )&
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_w_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_r,j-1,i), 28 )&
+                                     ), DIM = 1                                &
+                                        ) - 1
                 DO  k = nzb, nzt_topo_nestbc_r
 !
 !--                Wall for u on the south side, but not on the north side
-                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) )  .AND.        &
-                        ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )  THEN
+                   IF ( ( k_wall_u_ji > k_wall_u_ji_p )  .AND.                 &
+                        ( k_wall_u_ji == k_wall_u_ji_m ) )  THEN
                       inc        = 1
                       wall_index = j
                       CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1701,10 +2245,10 @@
 !
 !--                Wall for u on the north side, but not on the south side
-                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) )  .AND.        &
-                        ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
+                   IF ( ( k_wall_u_ji > k_wall_u_ji_m )  .AND.                  &
+                        ( k_wall_u_ji == k_wall_u_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  j+1
                       CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1717,10 +2261,10 @@
 !
 !--                Wall for w on the south side, but not on the north side
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j-1,i) ) )  THEN
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_p )  .AND.                 &
+                        ( k_wall_w_ji == k_wall_w_ji_m ) )  THEN
                       inc        =  1
                       wall_index =  j
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1733,10 +2277,10 @@
 !
 !--                Wall for w on the north side, but not on the south side
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_m )  .AND.                 &
+                        ( k_wall_w_ji == k_wall_w_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  j+1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
 
 !
@@ -1760,14 +2304,46 @@
 
              DO  i = nxl, nxr
+
+                k_wall_v_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,0,i), 27 ) &
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_v_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,0,i+1), 27 )&
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_v_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,0,i-1), 27 )&
+                                     ), DIM = 1                                &
+                                        ) - 1
+
+                k_wall_w_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,-1,i), 28 )&
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_w_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,-1,i+1), 28 )&
+                                 ), DIM = 1                                    &
+                                     ) - 1
+                k_wall_w_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_s,-1,i-1), 28 )&
+                                     ), DIM = 1                                &
+                                        ) - 1
                 DO  k = nzb, nzt_topo_nestbc_s
 !
 !--                Wall for v on the left side, but not on the right side
                    j  = 0
-                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.        &
-                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_p )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_m ) )  THEN
                       inc        =  1
                       wall_index =  i
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1781,10 +2357,10 @@
 !--                Wall for v on the right side, but not on the left side
                    j  = 0
-                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.        &
-                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_m )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  i+1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1798,10 +2374,10 @@
 !--                Wall for w on the left side, but not on the right side
                    j  = -1
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_p )  .AND.                 &
+                        ( k_wall_w_ji == k_wall_w_ji_m ) )  THEN
                       inc        =  1
                       wall_index =  i
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1815,10 +2391,10 @@
 !--                Wall for w on the right side, but not on the left side
                    j  = -1
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
+                   IF ( ( k_wall_w_ji > k_wall_w_ji_m )  .AND.                 &
+                        ( k_wall_w_ji == k_wall_w_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  i+1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1842,13 +2418,45 @@
 
              DO  i = nxl, nxr
+                k_wall_v_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i), 27 ) &
+                                 ), DIM = 1                                    &
+                                       ) - 1
+
+                k_wall_v_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i+1), 27 )&
+                                 ), DIM = 1                                    &
+                                       ) - 1
+                k_wall_v_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i-1), 27 )&
+                                     ), DIM = 1                                &
+                                       ) - 1
+
+                k_wall_w_ji    = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i), 28 )  &
+                                 ), DIM = 1                                    &
+                                       ) - 1
+                k_wall_w_ji_p  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i+1), 28 )&
+                                 ), DIM = 1                                    &
+                                       ) - 1
+                k_wall_w_ji_m  = MAXLOC(                                       &
+                            MERGE( 1, 0,                                       &
+                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_n,j,i-1), 28 )&
+                                     ), DIM = 1                                &
+                                       ) - 1
                 DO  k = nzb, nzt_topo_nestbc_n
 !
 !--                Wall for v on the left side, but not on the right side
-                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.        &
-                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_p )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_m ) )  THEN
                       inc        = 1
                       wall_index = i
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1861,10 +2469,10 @@
 !
 !--                Wall for v on the right side, but not on the left side
-                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.        &
-                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_m )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  i + 1
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1877,10 +2485,10 @@
 !
 !--                Wall for w on the left side, but not on the right side
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_p )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_m ) )  THEN
                       inc        = 1
                       wall_index = i
-                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1893,10 +2501,10 @@
 !
 !--                Wall for w on the right side, but not on the left side
-                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.        &
-                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
+                   IF ( ( k_wall_v_ji > k_wall_v_ji_m )  .AND.                 &
+                        ( k_wall_v_ji == k_wall_v_ji_p ) )  THEN
                       inc        = -1
                       wall_index =  i+1
                       CALL pmci_define_loglaw_correction_parameters( lc, lcr,   &
-                          k, i, inc, wall_index, z0(j,i), kb, direction, ncorr )
+                          k, i, inc, wall_index, z0_topo, kb, direction, ncorr )
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -2416,5 +3024,5 @@
 
 
-    SUBROUTINE pmci_init_tkefactor
+SUBROUTINE pmci_init_tkefactor
 
 !
@@ -2426,4 +3034,9 @@
 
        IMPLICIT NONE
+
+       INTEGER(iwp)        ::  k                     !: index variable along z
+       INTEGER(iwp)        ::  k_wall                !: topography-top index along z
+       INTEGER(iwp)        ::  kc                    !:
+
        REAL(wp), PARAMETER ::  cfw = 0.2_wp          !:
        REAL(wp), PARAMETER ::  c_tkef = 0.6_wp       !:
@@ -2434,8 +3047,5 @@
        REAL(wp)            ::  height                !:
        REAL(wp), PARAMETER ::  p13 = 1.0_wp/3.0_wp   !:
-       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !:
-       INTEGER(iwp)        ::  k                     !:
-       INTEGER(iwp)        ::  kc                    !:
-        
+       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !:        
 
        IF ( nest_bound_l )  THEN
@@ -2444,14 +3054,21 @@
           i = nxl - 1
           DO  j = nysg, nyng
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+
+             DO  k = k_wall + 1, nzt
+
                 kc     = kco(k+1)
                 glsf   = ( dx * dy * dzu(k) )**p13
                 glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
-                height = zu(k) - zu(nzb_s_inner(j,i))
+                height = zu(k) - zu(k_wall)
                 fw     = EXP( -cfw * height / glsf )
                 tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
                                               ( glsf / glsc )**p23 )
              ENDDO
-             tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0
+             tkefactor_l(k_wall,j) = c_tkef * fw0
           ENDDO
        ENDIF
@@ -2462,14 +3079,21 @@
           i = nxr + 1
           DO  j = nysg, nyng
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+
+             DO  k = k_wall + 1, nzt
+
                 kc     = kco(k+1)
                 glsf   = ( dx * dy * dzu(k) )**p13
                 glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                height = zu(k) - zu(nzb_s_inner(j,i))
+                height = zu(k) - zu(k_wall)
                 fw     = EXP( -cfw * height / glsf )
                 tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
                                               ( glsf / glsc )**p23 )
              ENDDO
-             tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0
+             tkefactor_r(k_wall,j) = c_tkef * fw0
           ENDDO
        ENDIF
@@ -2480,14 +3104,21 @@
           j = nys - 1
           DO  i = nxlg, nxrg
-             DO  k = nzb_s_inner(j,i) + 1, nzt
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+
+             DO  k = k_wall + 1, nzt
+
                 kc     = kco(k+1)
                 glsf   = ( dx * dy * dzu(k) )**p13
                 glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
-                height = zu(k) - zu(nzb_s_inner(j,i))
+                height = zu(k) - zu(k_wall)
                 fw     = EXP( -cfw*height / glsf )
                 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
                                               ( glsf / glsc )**p23 )
              ENDDO
-             tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0
+             tkefactor_s(k_wall,i) = c_tkef * fw0
           ENDDO
        ENDIF
@@ -2498,14 +3129,20 @@
           j = nyn + 1
           DO  i = nxlg, nxrg
-             DO  k = nzb_s_inner(j,i)+1, nzt
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall + 1, nzt
+
                 kc     = kco(k+1)
                 glsf   = ( dx * dy * dzu(k) )**p13
                 glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                height = zu(k) - zu(nzb_s_inner(j,i))
+                height = zu(k) - zu(k_wall)
                 fw     = EXP( -cfw * height / glsf )
-                tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
+                tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
                                               ( glsf / glsc )**p23 )
              ENDDO
-             tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0
+             tkefactor_n(k_wall,i) = c_tkef * fw0
           ENDDO
        ENDIF
@@ -2513,12 +3150,21 @@
        ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
        k = nzt
+
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
+!
+!--          Determine vertical index for local topography top
+             k_wall = MAXLOC(                                                  &
+                        MERGE( 1, 0,                                           &
+                               BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )      &
+                             ), DIM = 1                                        &
+                            ) - 1
+
              kc     = kco(k+1)
              glsf   = ( dx * dy * dzu(k) )**p13
              glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-             height = zu(k) - zu(nzb_s_inner(j,i))
+             height = zu(k) - zu(k_wall)
              fw     = EXP( -cfw * height / glsf )
-             tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *         &
+             tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
                                            ( glsf / glsc )**p23 )
           ENDDO
@@ -2783,4 +3429,5 @@
     INTEGER(iwp) ::  jcn        !:
     INTEGER(iwp) ::  jcs        !:
+    INTEGER(iwp) ::  k          !:
 
     REAL(wp) ::  waittime       !:
@@ -2804,15 +3451,15 @@
 !--    The interpolation.
        CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,    &
-                                   r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
+                                   r2yo, r1zo, r2zo, 'u' )
        CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,    &
-                                   r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
+                                   r2yv, r1zo, r2zo, 'v' )
        CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,    &
-                                   r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
+                                   r2yo, r1zw, r2zw, 'w' )
        CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,    &
-                                   r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
+                                   r2yo, r1zo, r2zo, 'e' )
 
        IF ( .NOT. neutral )  THEN
           CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,       &
-                                      r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
+                                      r1yo, r2yo, r1zo, r2zo, 's' )
        ENDIF
 
@@ -2820,19 +3467,15 @@
 
           CALL pmci_interp_tril_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo,   &
-                                      r2yo, r1zo, r2zo, nzb_s_inner, 's' )
+                                      r2yo, r1zo, r2zo, 's' )
 
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 !              CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,    &
-!                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner,   & 
-!                                          's' ) 
+!                                          r1yo, r2yo, r1zo, r2zo, 's' ) 
              CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,    &
-                                         r1yo, r2yo, r1zo, r2zo, nzb_s_inner,   & 
-                                         's' )
+                                         r1yo, r2yo, r1zo, r2zo, 's' )
 !             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,    &
-!                                         r1yo, r2yo, r1zo, r2zo, nzb_s_inner,   & 
-!                                         's' )   
+!                                         r1yo, r2yo, r1zo, r2zo, 's' )   
              CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,    &
-                                         r1yo, r2yo, r1zo, r2zo, nzb_s_inner,   & 
-                                         's' )
+                                         r1yo, r2yo, r1zo, r2zo, 's' )
           ENDIF
 
@@ -2841,5 +3484,5 @@
        IF ( passive_scalar )  THEN
           CALL pmci_interp_tril_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,   &
-                                      r2yo, r1zo, r2zo, nzb_s_inner, 's' )
+                                      r2yo, r1zo, r2zo, 's' )
        ENDIF
 
@@ -2851,12 +3494,22 @@
           DO   i = nxlg, nxrg
              DO   j = nysg, nyng
-                u(nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
-                v(nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
-                w(nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
-                e(nzb:nzb_s_inner(j,i),j,i)   = 0.0_wp
-                u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
-                v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
-                w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
-                e_p(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
+                DO  k = nzb, nzt
+                   u(k,j,i)   = MERGE( u(k,j,i), 0.0_wp,                       &
+                                       BTEST( wall_flags_0(k,j,i), 1 ) )
+                   v(k,j,i)   = MERGE( v(k,j,i), 0.0_wp,                       &
+                                       BTEST( wall_flags_0(k,j,i), 2 ) )
+                   w(k,j,i)   = MERGE( w(k,j,i), 0.0_wp,                       &
+                                       BTEST( wall_flags_0(k,j,i), 3 ) )
+                   e(k,j,i)   = MERGE( e(k,j,i), 0.0_wp,                       &
+                                       BTEST( wall_flags_0(k,j,i), 0 ) )
+                   u_p(k,j,i) = MERGE( u_p(k,j,i), 0.0_wp,                     &
+                                       BTEST( wall_flags_0(k,j,i), 1 ) )
+                   v_p(k,j,i) = MERGE( v_p(k,j,i), 0.0_wp,                     &
+                                       BTEST( wall_flags_0(k,j,i), 2 ) )
+                   w_p(k,j,i) = MERGE( w_p(k,j,i), 0.0_wp,                     &
+                                       BTEST( wall_flags_0(k,j,i), 3 ) )
+                   e_p(k,j,i) = MERGE( e_p(k,j,i), 0.0_wp,                     &
+                                       BTEST( wall_flags_0(k,j,i), 0 ) )
+                ENDDO
              ENDDO
           ENDDO
@@ -2869,5 +3522,5 @@
 
     SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,     &
-                                     r1z, r2z, kb, var )
+                                     r1z, r2z, var )
 !
 !--    Interpolation of the internal values for the child-domain initialization
@@ -2881,18 +3534,20 @@
        INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !:
        INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
-       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !:
-
-       INTEGER(iwp) ::  i      !:
-       INTEGER(iwp) ::  ib     !:
-       INTEGER(iwp) ::  ie     !:
-       INTEGER(iwp) ::  j      !:
-       INTEGER(iwp) ::  jb     !:
-       INTEGER(iwp) ::  je     !:
-       INTEGER(iwp) ::  k      !:
-       INTEGER(iwp) ::  k1     !:
-       INTEGER(iwp) ::  kbc    !:
-       INTEGER(iwp) ::  l      !:
-       INTEGER(iwp) ::  m      !:
-       INTEGER(iwp) ::  n      !:
+
+       INTEGER(iwp) ::  flag_nr  !: Number of flag array to mask topography on respective u/v/w or s grid
+       INTEGER(iwp) ::  flag_nr2 !: Number of flag array to indicate vertical index of topography top on respective u/v/w or s grid
+       INTEGER(iwp) ::  i        !:
+       INTEGER(iwp) ::  ib       !:
+       INTEGER(iwp) ::  ie       !:
+       INTEGER(iwp) ::  j        !:
+       INTEGER(iwp) ::  jb       !:
+       INTEGER(iwp) ::  je       !:
+       INTEGER(iwp) ::  k        !:
+       INTEGER(iwp) ::  k_wall   !:
+       INTEGER(iwp) ::  k1       !:
+       INTEGER(iwp) ::  kbc      !:
+       INTEGER(iwp) ::  l        !:
+       INTEGER(iwp) ::  m        !:
+       INTEGER(iwp) ::  n        !:
 
        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
@@ -2914,4 +3569,5 @@
        REAL(wp) ::  logzuc1    !:
        REAL(wp) ::  zuc1       !:
+       REAL(wp) ::  z0_topo    !:  roughness at vertical walls
 
 
@@ -2945,8 +3601,23 @@
        ENDIF
 !
+!--    Determine number of flag array to be used to mask topography
+       IF ( var == 'u' )  THEN
+          flag_nr  = 1
+          flag_nr2 = 14
+       ELSEIF ( var == 'v' )  THEN
+          flag_nr  = 2
+          flag_nr2 = 16
+       ELSEIF ( var == 'w' )  THEN
+          flag_nr  = 3
+          flag_nr2 = 18
+       ELSE
+          flag_nr  = 0
+          flag_nr2 = 12
+       ENDIF
+!
 !--    Trilinear interpolation.
        DO  i = ib, ie
           DO  j = jb, je
-             DO  k = kb(j,i), nzt + 1
+             DO  k = nzb, nzt + 1
                 l = ic(i)
                 m = jc(j)
@@ -2970,27 +3641,35 @@
 !--    too.
        IF ( var == 'u' .OR. var == 'v' )  THEN
+          z0_topo = roughness_length
           DO  i = ib, nxr
              DO  j = jb, nyn
+!
+!--             Determine vertical index of topography top at grid point (j,i)
+                k_wall = MAXLOC(                                               &
+                      MERGE( 1, 0,                                             &
+                             BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )  &
+                           ), DIM = 1                                          &
+                               ) - 1
+!
+!--             kbc is the first coarse-grid point above the surface
                 kbc = 1
-!
-!--             kbc is the first coarse-grid point above the surface
-                DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )
+                DO  WHILE ( cg%zu(kbc) < zu(k_wall) )
                    kbc = kbc + 1
                 ENDDO
                 zuc1 = cg%zu(kbc)
-                k1   = kb(j,i) + 1
+                k1   = k_wall + 1
                 DO  WHILE ( zu(k1) < zuc1 )
                    k1 = k1 + 1
                 ENDDO
-                logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
-
-                k = kb(j,i) + 1
+                logzuc1 = LOG( ( zu(k1) - zu(k_wall) ) / z0_topo )
+
+                k = k_wall + 1
                 DO  WHILE ( zu(k) < zuc1 )
-                   logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) /     &
+                   logratio = ( LOG( ( zu(k) - zu(k_wall) ) / z0_topo ) ) /     &
                                 logzuc1
                    f(k,j,i) = logratio * f(k1,j,i)
                    k  = k + 1
                 ENDDO
-                f(kb(j,i),j,i) = 0.0_wp
+                f(k_wall,j,i) = 0.0_wp
              ENDDO
           ENDDO
@@ -3000,7 +3679,15 @@
           DO  i = ib, nxr
               DO  j = jb, nyn
-                f(kb(j,i),j,i) = 0.0_wp
-             ENDDO
-          ENDDO
+!
+!--              Determine vertical index of topography top at grid point (j,i)
+                 k_wall = MAXLOC(                                              &
+                      MERGE( 1, 0,                                             &
+                             BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )  &
+                           ), DIM = 1                                          &
+                                ) - 1
+
+                 f(k_wall,j,i) = 0.0_wp
+              ENDDO
+           ENDDO
 
        ENDIF
@@ -3137,6 +3824,8 @@
        innor = dy
        DO   j = nys, nyn
-          DO   k = nzb_u_inner(j,i)+1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
+          DO   k = nzb+1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
           ENDDO
        ENDDO
@@ -3147,6 +3836,8 @@
        innor = -dy
        DO   j = nys, nyn
-          DO   k = nzb_u_inner(j,i)+1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
+          DO   k = nzb+1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 1 ) )
           ENDDO
        ENDDO
@@ -3170,6 +3861,8 @@
        innor = dx
        DO   i = nxl, nxr
-          DO   k = nzb_v_inner(j,i)+1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
+          DO   k = nzb+1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
           ENDDO
        ENDDO
@@ -3180,6 +3873,8 @@
        innor = -dx
        DO   i = nxl, nxr
-          DO   k = nzb_v_inner(j,i)+1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
+          DO   k = nzb+1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
+                                     * MERGE( 1.0_wp, 0.0_wp,                  &
+                                              BTEST( wall_flags_0(k,j,i), 2 ) )
           ENDDO
        ENDDO
@@ -3340,6 +4035,7 @@
     INTEGER(iwp) ::  child_id    !:
     INTEGER(iwp) ::  i           !:
+    INTEGER(iwp) ::  ierr        !:
     INTEGER(iwp) ::  j           !:
-    INTEGER(iwp) ::  ierr        !:
+    INTEGER(iwp) ::  k           !:
     INTEGER(iwp) ::  m           !:
 
@@ -3374,8 +4070,13 @@
              DO   i = nxlg, nxrg
                 DO   j = nysg, nyng
-                   u(nzb:nzb_u_inner(j,i),j,i)  = 0.0_wp
-                   v(nzb:nzb_v_inner(j,i),j,i)  = 0.0_wp
-                   w(nzb:nzb_w_inner(j,i),j,i)  = 0.0_wp
-                   e(nzb:nzb_s_inner(j,i),j,i)  = 0.0_wp
+                   DO  k = nzb, nzt+1
+                      u(k,j,i)  = MERGE( u(k,j,i), 0.0_wp,                     &
+                                         BTEST( wall_flags_0(k,j,i), 1 ) )
+                      v(k,j,i)  = MERGE( v(k,j,i), 0.0_wp,                     &
+                                         BTEST( wall_flags_0(k,j,i), 2 ) )
+                      w(k,j,i)  = MERGE( w(k,j,i), 0.0_wp,                     &
+                                         BTEST( wall_flags_0(k,j,i), 3 ) )
+                      e(k,j,i)  = MERGE( e(k,j,i), 0.0_wp,                     &
+                                         BTEST( wall_flags_0(k,j,i), 0 ) )
 !
 !--                TO_DO: zero setting of temperature within topography creates
@@ -3385,4 +4086,5 @@
 !                      q(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
 !                   ENDIF
+                   ENDDO
                 ENDDO
              ENDDO
@@ -3463,20 +4165,20 @@
           IF ( nest_bound_l )  THEN
              CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_u_inner,     &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
                                        logc_u_l, logc_ratio_u_l,                &
                                        nzt_topo_nestbc_l, 'l', 'u' )
 
              CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
-                                       r1yv, r2yv, r1zo, r2zo, nzb_v_inner,     &
+                                       r1yv, r2yv, r1zo, r2zo,                  &
                                        logc_v_l, logc_ratio_v_l,                &
                                        nzt_topo_nestbc_l, 'l', 'v' )
 
              CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
-                                       r1yo, r2yo, r1zw, r2zw, nzb_w_inner,     &
+                                       r1yo, r2yo, r1zw, r2zw,                  &
                                        logc_w_l, logc_ratio_w_l,                &
                                        nzt_topo_nestbc_l, 'l', 'w' )
 
              CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_s_inner,     &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
                                        logc_u_l, logc_ratio_u_l,                &
                                        nzt_topo_nestbc_l, 'l', 'e' )
@@ -3484,5 +4186,5 @@
              IF ( .NOT. neutral )  THEN
                 CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,   &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner,  &
+                                          r1yo, r2yo, r1zo, r2zo,               &
                                           logc_u_l, logc_ratio_u_l,             &
                                           nzt_topo_nestbc_l, 'l', 's' )
@@ -3492,5 +4194,5 @@
 
                 CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,    &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner,  &
+                                          r1yo, r2yo, r1zo, r2zo,               &
                                           logc_u_l, logc_ratio_u_l,             &
                                           nzt_topo_nestbc_l, 'l', 's' )
@@ -3500,5 +4202,5 @@
 !                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo, r2xo,&
 !                                              r1yo, r2yo, r1zo, r2zo,            &
-!                                              nzb_s_inner, logc_u_l,             &
+!                                              logc_u_l,                          &
 !                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
 !                                              'l', 's' )  
@@ -3506,5 +4208,5 @@
                    CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo, r2xo,&
                                              r1yo, r2yo, r1zo, r2zo,            &
-                                             nzb_s_inner, logc_u_l,             &
+                                             logc_u_l,                          &
                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
                                              'l', 's' ) 
@@ -3512,5 +4214,5 @@
 !                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, r2xo,&
 !                                              r1yo, r2yo, r1zo, r2zo,            &
-!                                              nzb_s_inner, logc_u_l,             &
+!                                              logc_u_l,                          &
 !                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
 !                                              'l', 's' ) 
@@ -3518,5 +4220,5 @@
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo, r2xo,&
                                              r1yo, r2yo, r1zo, r2zo,            &
-                                             nzb_s_inner, logc_u_l,             &
+                                             logc_u_l,                          &
                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
                                              'l', 's' )             
@@ -3527,5 +4229,5 @@
              IF ( passive_scalar )  THEN
                 CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,     &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner,  &
+                                          r1yo, r2yo, r1zo, r2zo,               &
                                           logc_u_l, logc_ratio_u_l,             &
                                           nzt_topo_nestbc_l, 'l', 's' )
@@ -3533,24 +4235,22 @@
 
              IF ( TRIM( nesting_mode ) == 'one-way' )  THEN
-
-                CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
-                CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
-                CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
-                CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
+                CALL pmci_extrap_ifoutflow_lr( u, 'l', 'u' )
+                CALL pmci_extrap_ifoutflow_lr( v, 'l', 'v' )
+                CALL pmci_extrap_ifoutflow_lr( w, 'l', 'w' )
+                CALL pmci_extrap_ifoutflow_lr( e, 'l', 'e' )
 
                 IF ( .NOT. neutral )  THEN
-                   CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( pt, 'l', 's' )
                 ENDIF
 
                 IF ( humidity )  THEN
-
-                   CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( q, 'l', 's' )
 
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
-!                       CALL pmci_extrap_ifoutflow_lr( qc, nzb_s_inner, 'l', 's' )
-                      CALL pmci_extrap_ifoutflow_lr( qr, nzb_s_inner, 'l', 's' )
-!                      CALL pmci_extrap_ifoutflow_lr( nc, nzb_s_inner, 'l', 's' )
-                      CALL pmci_extrap_ifoutflow_lr( nr, nzb_s_inner, 'l', 's' )
+!                       CALL pmci_extrap_ifoutflow_lr( qc, 'l', 's' )
+                      CALL pmci_extrap_ifoutflow_lr( qr, 'l', 's' )
+!                      CALL pmci_extrap_ifoutflow_lr( nc, 'l', 's' )
+                      CALL pmci_extrap_ifoutflow_lr( nr, 'l', 's' )
 
                    ENDIF
@@ -3559,5 +4259,5 @@
 
                 IF ( passive_scalar )  THEN
-                   CALL pmci_extrap_ifoutflow_lr( s, nzb_s_inner, 'l', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( s, 'l', 's' )
                 ENDIF
 
@@ -3569,38 +4269,39 @@
    !--    Right border pe
           IF ( nest_bound_r )  THEN
-
-             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_u_inner,    &
-                                       logc_u_r, logc_ratio_u_r,               &
+             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_r, logc_ratio_u_r,                &
                                        nzt_topo_nestbc_r, 'r', 'u' )
 
-             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
-                                       r1yv, r2yv, r1zo, r2zo, nzb_v_inner,    &
-                                       logc_v_r, logc_ratio_v_r,               &
+             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
+                                       r1yv, r2yv, r1zo, r2zo,                  &
+                                       logc_v_r, logc_ratio_v_r,                &
                                        nzt_topo_nestbc_r, 'r', 'v' )
 
-             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zw, r2zw, nzb_w_inner,    &
-                                       logc_w_r, logc_ratio_w_r,               &
+             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zw, r2zw,                  &
+                                       logc_w_r, logc_ratio_w_r,                &
                                        nzt_topo_nestbc_r, 'r', 'w' )
 
-             CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,     &
-                                       r1yo,r2yo, r1zo, r2zo, nzb_s_inner,     &
-                                       logc_u_r, logc_ratio_u_r,               &
+             CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
+                                       r1yo,r2yo, r1zo, r2zo,                   &
+                                       logc_u_r, logc_ratio_u_r,                &
                                        nzt_topo_nestbc_r, 'r', 'e' )
 
+
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_r, logc_ratio_u_r,            &
+                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,   &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_r, logc_ratio_u_r,             &
                                           nzt_topo_nestbc_r, 'r', 's' )
+
              ENDIF
 
              IF ( humidity )  THEN
-
-                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_r, logc_ratio_u_r,            &
+                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,    &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_r, logc_ratio_u_r,             &
                                           nzt_topo_nestbc_r, 'r', 's' )
+
 
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -3608,5 +4309,5 @@
 !                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              nzb_s_inner, logc_u_r,            &
+!                                              logc_u_r,                         &
 !                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
 !                                              'r', 's' ) 
@@ -3614,5 +4315,5 @@
                    CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-                                             nzb_s_inner, logc_u_r,            &
+                                             logc_u_r,                         &
                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
                                              'r', 's' ) 
@@ -3620,5 +4321,5 @@
 !                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              nzb_s_inner, logc_u_r,            &
+!                                              logc_u_r,                         &
 !                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
 !                                              'r', 's' ) 
@@ -3626,5 +4327,5 @@
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-                                             nzb_s_inner, logc_u_r,            &
+                                             logc_u_r,                         &
                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
                                              'r', 's' ) 
@@ -3636,29 +4337,28 @@
              IF ( passive_scalar )  THEN
                 CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
+                                          r1yo, r2yo, r1zo, r2zo,              &
                                           logc_u_r, logc_ratio_u_r,            &
                                           nzt_topo_nestbc_r, 'r', 's' )
+
              ENDIF
 
              IF ( TRIM( nesting_mode ) == 'one-way' )  THEN
-
-                CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
-                CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'r', 'v' )
-                CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'r', 'w' )
-                CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'r', 'e' )
+                CALL pmci_extrap_ifoutflow_lr( u, 'r', 'u' )
+                CALL pmci_extrap_ifoutflow_lr( v, 'r', 'v' )
+                CALL pmci_extrap_ifoutflow_lr( w, 'r', 'w' )
+                CALL pmci_extrap_ifoutflow_lr( e, 'r', 'e' )
 
                 IF ( .NOT. neutral )  THEN
-                   CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'r', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( pt, 'r', 's' )
                 ENDIF
 
                 IF ( humidity )  THEN
-
-                   CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'r', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( q, 'r', 's' )
 
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_lr( qc, nzb_s_inner, 'r', 's' )
-                      CALL pmci_extrap_ifoutflow_lr( qr, nzb_s_inner, 'r', 's' )
-!                      CALL pmci_extrap_ifoutflow_lr( nc, nzb_s_inner, 'r', 's' )  
-                      CALL pmci_extrap_ifoutflow_lr( nr, nzb_s_inner, 'r', 's' )
+!                       CALL pmci_extrap_ifoutflow_lr( qc, 'r', 's' )
+                      CALL pmci_extrap_ifoutflow_lr( qr, 'r', 's' )
+!                      CALL pmci_extrap_ifoutflow_lr( nc, 'r', 's' )  
+                      CALL pmci_extrap_ifoutflow_lr( nr, 'r', 's' )
                    ENDIF
 
@@ -3666,5 +4366,5 @@
 
                 IF ( passive_scalar )  THEN
-                   CALL pmci_extrap_ifoutflow_lr( s, nzb_s_inner, 'r', 's' )
+                   CALL pmci_extrap_ifoutflow_lr( s, 'r', 's' )
                 ENDIF
              ENDIF
@@ -3675,36 +4375,32 @@
    !--    South border pe
           IF ( nest_bound_s )  THEN
-             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_u_inner,    &
-                                       logc_u_s, logc_ratio_u_s,               &
+             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_s, logc_ratio_u_s,                &
                                        nzt_topo_nestbc_s, 's', 'u' )
-
-             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
-                                       r1yv, r2yv, r1zo, r2zo, nzb_v_inner,    &
-                                       logc_v_s, logc_ratio_v_s,               &
+             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
+                                       r1yv, r2yv, r1zo, r2zo,                  &
+                                       logc_v_s, logc_ratio_v_s,                &
                                        nzt_topo_nestbc_s, 's', 'v' )
-
-             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zw, r2zw, nzb_w_inner,    &
-                                       logc_w_s, logc_ratio_w_s,               &
+             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zw, r2zw,                  &
+                                       logc_w_s, logc_ratio_w_s,                &
                                        nzt_topo_nestbc_s, 's','w' )
-
-             CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_s_inner,    &
-                                       logc_u_s, logc_ratio_u_s,               &
+             CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_s, logc_ratio_u_s,                &
                                        nzt_topo_nestbc_s, 's', 'e' )
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_s, logc_ratio_u_s,            &
+                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,   &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_s, logc_ratio_u_s,             &
                                           nzt_topo_nestbc_s, 's', 's' )
              ENDIF
 
              IF ( humidity )  THEN
-
-                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
-                                          r1yo,r2yo, r1zo, r2zo, nzb_s_inner,  &
-                                          logc_u_s, logc_ratio_u_s,            &
+                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,    &
+                                          r1yo,r2yo, r1zo, r2zo,                &
+                                          logc_u_s, logc_ratio_u_s,             &
                                           nzt_topo_nestbc_s, 's', 's' )
 
@@ -3713,5 +4409,5 @@
 !                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-!                                              nzb_s_inner, logc_u_s,            &
+!                                              logc_u_s,                         &
 !                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
 !                                              's', 's' )
@@ -3719,5 +4415,5 @@
                    CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-                                             nzb_s_inner, logc_u_s,            &
+                                             logc_u_s,                         &
                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
                                              's', 's' )
@@ -3725,5 +4421,5 @@
 !                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-!                                              nzb_s_inner, logc_u_s,            &
+!                                              logc_u_s,                         &
 !                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
 !                                              's', 's' )
@@ -3731,5 +4427,5 @@
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-                                             nzb_s_inner, logc_u_s,            &
+                                             logc_u_s,                         &
                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
                                              's', 's' )
@@ -3740,30 +4436,28 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
-                                          r1yo,r2yo, r1zo, r2zo, nzb_s_inner,  &
-                                          logc_u_s, logc_ratio_u_s,            &
+                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,     &
+                                          r1yo,r2yo, r1zo, r2zo,                &
+                                          logc_u_s, logc_ratio_u_s,             &
                                           nzt_topo_nestbc_s, 's', 's' )
              ENDIF
 
              IF ( TRIM( nesting_mode ) == 'one-way' )  THEN
-
-                CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 's', 'u' )
-                CALL pmci_extrap_ifoutflow_sn( v, nzb_v_inner, 's', 'v' )
-                CALL pmci_extrap_ifoutflow_sn( w, nzb_w_inner, 's', 'w' )
-                CALL pmci_extrap_ifoutflow_sn( e, nzb_s_inner, 's', 'e' )
+                CALL pmci_extrap_ifoutflow_sn( u, 's', 'u' )
+                CALL pmci_extrap_ifoutflow_sn( v, 's', 'v' )
+                CALL pmci_extrap_ifoutflow_sn( w, 's', 'w' )
+                CALL pmci_extrap_ifoutflow_sn( e, 's', 'e' )
 
                 IF ( .NOT. neutral )  THEN
-                   CALL pmci_extrap_ifoutflow_sn( pt,nzb_s_inner, 's', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( pt, 's', 's' )
                 ENDIF
 
                 IF ( humidity )  THEN
-
-                   CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 's', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( q,  's', 's' )
 
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_sn( qc, nzb_s_inner, 's', 's' ) 
-                      CALL pmci_extrap_ifoutflow_sn( qr, nzb_s_inner, 's', 's' )      
-!                       CALL pmci_extrap_ifoutflow_sn( nc, nzb_s_inner, 's', 's' ) 
-                      CALL pmci_extrap_ifoutflow_sn( nr, nzb_s_inner, 's', 's' ) 
+!                       CALL pmci_extrap_ifoutflow_sn( qc, 's', 's' ) 
+                      CALL pmci_extrap_ifoutflow_sn( qr, 's', 's' )      
+!                       CALL pmci_extrap_ifoutflow_sn( nc, 's', 's' ) 
+                      CALL pmci_extrap_ifoutflow_sn( nr, 's', 's' ) 
 
                    ENDIF
@@ -3772,5 +4466,5 @@
 
                 IF ( passive_scalar )  THEN
-                   CALL pmci_extrap_ifoutflow_sn( s, nzb_s_inner, 's', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( s,  's', 's' )
                 ENDIF
 
@@ -3782,34 +4476,35 @@
    !--    North border pe
           IF ( nest_bound_n )  THEN
-
-             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_u_inner,    &
-                                       logc_u_n, logc_ratio_u_n,               &
+             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_n, logc_ratio_u_n,                &
                                        nzt_topo_nestbc_n, 'n', 'u' )
-             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
-                                       r1yv, r2yv, r1zo, r2zo, nzb_v_inner,    &
-                                       logc_v_n, logc_ratio_v_n,               & 
+
+             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
+                                       r1yv, r2yv, r1zo, r2zo,                  &
+                                       logc_v_n, logc_ratio_v_n,                & 
                                        nzt_topo_nestbc_n, 'n', 'v' )
-             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zw, r2zw, nzb_w_inner,    &
-                                       logc_w_n, logc_ratio_w_n,               &
+
+             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zw, r2zw,                  &
+                                       logc_w_n, logc_ratio_w_n,                &
                                        nzt_topo_nestbc_n, 'n', 'w' )
-             CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zo, r2zo, nzb_s_inner,    &
-                                       logc_u_n, logc_ratio_u_n,               &
+
+             CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_n, logc_ratio_u_n,                &
                                        nzt_topo_nestbc_n, 'n', 'e' )
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_n, logc_ratio_u_n,            &
+                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,   &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_n, logc_ratio_u_n,             &
                                           nzt_topo_nestbc_n, 'n', 's' )
              ENDIF
 
              IF ( humidity )  THEN
-
-                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_n, logc_ratio_u_n,            &
+                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,    &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_n, logc_ratio_u_n,             &
                                           nzt_topo_nestbc_n, 'n', 's' )
 
@@ -3818,5 +4513,5 @@
 !                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              nzb_s_inner, logc_u_n,            &
+!                                              logc_u_n,                         &
 !                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
 !                                              'n', 's' )
@@ -3824,5 +4519,5 @@
                    CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-                                             nzb_s_inner, logc_u_n,            &
+                                             logc_u_n,                         &
                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
                                              'n', 's' )
@@ -3830,5 +4525,5 @@
 !                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              nzb_s_inner, logc_u_n,            &
+!                                              logc_u_n,                         &
 !                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
 !                                              'n', 's' )
@@ -3836,5 +4531,5 @@
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-                                             nzb_s_inner, logc_u_n,            &
+                                             logc_u_n,                         &
                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
                                              'n', 's' )
@@ -3845,30 +4540,28 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
-                                          r1yo, r2yo, r1zo, r2zo, nzb_s_inner, &
-                                          logc_u_n, logc_ratio_u_n,            &
+                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,     &
+                                          r1yo, r2yo, r1zo, r2zo,               &
+                                          logc_u_n, logc_ratio_u_n,             &
                                           nzt_topo_nestbc_n, 'n', 's' )
              ENDIF
 
              IF ( TRIM( nesting_mode ) == 'one-way' )  THEN
-
-                CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 'n', 'u' )
-                CALL pmci_extrap_ifoutflow_sn( v, nzb_v_inner, 'n', 'v' )
-                CALL pmci_extrap_ifoutflow_sn( w, nzb_w_inner, 'n', 'w' )
-                CALL pmci_extrap_ifoutflow_sn( e, nzb_s_inner, 'n', 'e' )
+                CALL pmci_extrap_ifoutflow_sn( u, 'n', 'u' )
+                CALL pmci_extrap_ifoutflow_sn( v, 'n', 'v' )
+                CALL pmci_extrap_ifoutflow_sn( w, 'n', 'w' )
+                CALL pmci_extrap_ifoutflow_sn( e, 'n', 'e' )
 
                 IF ( .NOT. neutral )  THEN
-                   CALL pmci_extrap_ifoutflow_sn( pt,nzb_s_inner, 'n', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( pt, 'n', 's' )
                 ENDIF
 
                 IF ( humidity )  THEN
-
-                   CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 'n', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( q,  'n', 's' )
 
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_sn( qc, nzb_s_inner, 'n', 's' )
-                      CALL pmci_extrap_ifoutflow_sn( qr, nzb_s_inner, 'n', 's' )
-!                       CALL pmci_extrap_ifoutflow_sn( nc, nzb_s_inner, 'n', 's' )
-                      CALL pmci_extrap_ifoutflow_sn( nr, nzb_s_inner, 'n', 's' )
+!                       CALL pmci_extrap_ifoutflow_sn( qc, 'n', 's' )
+                      CALL pmci_extrap_ifoutflow_sn( qr, 'n', 's' )
+!                       CALL pmci_extrap_ifoutflow_sn( nc, 'n', 's' )
+                      CALL pmci_extrap_ifoutflow_sn( nr, 'n', 's' )
                    ENDIF
 
@@ -3876,5 +4569,5 @@
 
                 IF ( passive_scalar )  THEN
-                   CALL pmci_extrap_ifoutflow_sn( s, nzb_s_inner, 'n', 's' )
+                   CALL pmci_extrap_ifoutflow_sn( s,  'n', 's' )
                 ENDIF
 
@@ -3883,5 +4576,5 @@
           ENDIF
 
-       ENDIF       !: IF ( nesting_mode /= 'vertical' )
+       ENDIF       ! IF ( nesting_mode /= 'vertical' )
 
 !
@@ -4016,5 +4709,5 @@
 
    SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z,  &
-                                   r2z, kb, logc, logc_ratio, nzt_topo_nestbc,  &
+                                   r2z, logc, logc_ratio, nzt_topo_nestbc,      &
                                    edge, var )
 !
@@ -4040,5 +4733,4 @@
       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic     !:
       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc     !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb     !:
       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !:
       INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nys:nyn),                 &
@@ -4049,22 +4741,25 @@
       CHARACTER(LEN=1), INTENT(IN) ::  var    !:
 
-      INTEGER(iwp) ::  i       !:
-      INTEGER(iwp) ::  ib      !:
-      INTEGER(iwp) ::  ibgp    !:
-      INTEGER(iwp) ::  iw      !:
-      INTEGER(iwp) ::  j       !:
-      INTEGER(iwp) ::  jco     !:
-      INTEGER(iwp) ::  jcorr   !:
-      INTEGER(iwp) ::  jinc    !:
-      INTEGER(iwp) ::  jw      !:
-      INTEGER(iwp) ::  j1      !:
-      INTEGER(iwp) ::  k       !:
-      INTEGER(iwp) ::  kco     !:
-      INTEGER(iwp) ::  kcorr   !:
-      INTEGER(iwp) ::  k1      !:
-      INTEGER(iwp) ::  l       !:
-      INTEGER(iwp) ::  m       !:
-      INTEGER(iwp) ::  n       !:
-      INTEGER(iwp) ::  kbc     !:
+      INTEGER(iwp) ::  flag_nr  !: Number of flag array to mask topography on respective u/v/w or s grid
+      INTEGER(iwp) ::  flag_nr2 !: Number of flag array to indicate vertical index of topography top on respective u/v/w or s grid
+      INTEGER(iwp) ::  i        !:
+      INTEGER(iwp) ::  ib       !:
+      INTEGER(iwp) ::  ibgp     !:
+      INTEGER(iwp) ::  iw       !:
+      INTEGER(iwp) ::  j        !:
+      INTEGER(iwp) ::  jco      !:
+      INTEGER(iwp) ::  jcorr    !:
+      INTEGER(iwp) ::  jinc     !:
+      INTEGER(iwp) ::  jw       !:
+      INTEGER(iwp) ::  j1       !:
+      INTEGER(iwp) ::  k        !:
+      INTEGER(iwp) ::  k_wall   !: vertical index of topography top
+      INTEGER(iwp) ::  kco      !:
+      INTEGER(iwp) ::  kcorr    !:
+      INTEGER(iwp) ::  k1       !:
+      INTEGER(iwp) ::  l        !:
+      INTEGER(iwp) ::  m        !:
+      INTEGER(iwp) ::  n        !:
+      INTEGER(iwp) ::  kbc      !:
       
       REAL(wp) ::  coarse_dx   !:
@@ -4094,7 +4789,22 @@
          ib = nxr + 2
       ENDIF
+!
+!--    Determine number of flag array to be used to mask topography
+       IF ( var == 'u' )  THEN
+          flag_nr  = 1
+          flag_nr2 = 14
+       ELSEIF ( var == 'v' )  THEN
+          flag_nr  = 2
+          flag_nr2 = 16
+       ELSEIF ( var == 'w' )  THEN
+          flag_nr  = 3
+          flag_nr2 = 18
+       ELSE
+          flag_nr  = 0
+          flag_nr2 = 12
+       ENDIF
       
       DO  j = nys, nyn+1
-         DO  k = kb(j,i), nzt+1
+         DO  k = nzb, nzt+1
             l = ic(i)
             m = jc(j)
@@ -4119,5 +4829,13 @@
       IF ( var == 'u' .OR. var == 'v' )  THEN           
          DO  j = nys, nyn
-            k = kb(j,i)+1
+!
+!--         Determine vertical index of topography top at grid point (j,i)
+            k_wall = MAXLOC(                                                   &
+                        MERGE( 1, 0,                                           &
+                               BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )&
+                             ), DIM = 1                                        &
+                           ) - 1
+
+            k = k_wall+1
             IF ( ( logc(1,k,j) /= 0 )  .AND.  ( logc(2,k,j) == 0 ) )  THEN
                k1 = logc(1,k,j)
@@ -4141,5 +4859,12 @@
 !--         Solid surface only on south/north side of the node                   
             DO  j = nys, nyn
-               DO  k = kb(j,i)+1, nzt_topo_nestbc
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall+1, nzt_topo_nestbc
                   IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) == 0 ) )  THEN
 !
@@ -4162,5 +4887,12 @@
          IF ( var == 'u' )  THEN
             DO  j = nys, nyn
-               k = kb(j,i) + 1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               k = k_wall + 1
                IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) /= 0 ) )  THEN
                   k1   = logc(1,k,j)                 
@@ -4190,5 +4922,12 @@
          IF ( edge == 'l' )  THEN
             DO  j = nys, nyn + 1
-               DO  k = kb(j,i), nzt + 1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall, nzt + 1
                   f(k,j,i) = tkefactor_l(k,j) * f(k,j,i)
                ENDDO
@@ -4196,5 +4935,12 @@
          ELSEIF ( edge == 'r' )  THEN           
             DO  j = nys, nyn+1
-               DO  k = kb(j,i), nzt+1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
                ENDDO
@@ -4220,5 +4966,5 @@
 
    SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z,  &
-                                   r2z, kb, logc, logc_ratio,                   &
+                                   r2z, logc, logc_ratio,                   &
                                    nzt_topo_nestbc, edge, var )
 
@@ -4245,5 +4991,4 @@
       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !:
       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !:
       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
       INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nxl:nxr),                 &
@@ -4254,4 +4999,6 @@
       CHARACTER(LEN=1), INTENT(IN) ::  var    !:
       
+      INTEGER(iwp) ::  flag_nr  !: Number of flag array to mask topography on respective u/v/w or s grid
+      INTEGER(iwp) ::  flag_nr2 !: Number of flag array to indicate vertical index of topography top on respective u/v/w or s grid
       INTEGER(iwp) ::  i       !:
       INTEGER(iwp) ::  iinc    !:
@@ -4263,4 +5010,5 @@
       INTEGER(iwp) ::  jbgp    !:
       INTEGER(iwp) ::  k       !:
+      INTEGER(iwp) ::  k_wall   !: vertical index of topography top
       INTEGER(iwp) ::  kcorr   !:
       INTEGER(iwp) ::  kco     !:
@@ -4297,6 +5045,29 @@
       ENDIF
 
+!
+!--    Determine number of flag array to be used to mask topography
+       IF ( var == 'u' )  THEN
+          flag_nr  = 1
+          flag_nr2 = 14
+       ELSEIF ( var == 'v' )  THEN
+          flag_nr  = 2
+          flag_nr2 = 16
+       ELSEIF ( var == 'w' )  THEN
+          flag_nr  = 3
+          flag_nr2 = 18
+       ELSE
+          flag_nr  = 0
+          flag_nr2 = 12
+       ENDIF
+
       DO  i = nxl, nxr+1
-         DO  k = kb(j,i), nzt+1
+!
+!--      Determine vertical index of topography top at grid point (j,i)
+         k_wall = MAXLOC(                                                      &
+                        MERGE( 1, 0,                                           &
+                               BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )&
+                             ), DIM = 1                                        &
+                           ) - 1
+         DO  k = k_wall, nzt+1
             l = ic(i)
             m = jc(j)
@@ -4321,5 +5092,13 @@
       IF ( var == 'u'  .OR.  var == 'v' )  THEN           
          DO  i = nxl, nxr
-            k = kb(j,i) + 1
+!
+!--         Determine vertical index of topography top at grid point (j,i)
+            k_wall = MAXLOC(                                                   &
+                        MERGE( 1, 0,                                           &
+                               BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )&
+                             ), DIM = 1                                        &
+                           ) - 1
+
+            k = k_wall + 1
             IF ( ( logc(1,k,i) /= 0 )  .AND.  ( logc(2,k,i) == 0 ) )  THEN
                k1 = logc(1,k,i)
@@ -4342,5 +5121,12 @@
          IF ( var == 'v' .OR. var == 'w' )  THEN
             DO  i = nxl, nxr
-               DO  k = kb(j,i), nzt_topo_nestbc
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall, nzt_topo_nestbc
 !
 !--               Solid surface only on left/right side of the node           
@@ -4365,5 +5151,12 @@
          IF ( var == 'v' )  THEN
             DO  i = nxl, nxr
-               k = kb(j,i) + 1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               k = k_wall + 1
                IF ( ( logc(2,k,i) /= 0 )  .AND.  ( logc(1,k,i) /= 0 ) )  THEN
                   k1   = logc(1,k,i)         
@@ -4393,5 +5186,12 @@
          IF ( edge == 's' )  THEN
             DO  i = nxl, nxr + 1
-               DO  k = kb(j,i), nzt+1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
                ENDDO
@@ -4399,5 +5199,12 @@
          ELSEIF ( edge == 'n' )  THEN
             DO  i = nxl, nxr + 1
-               DO  k = kb(j,i), nzt+1
+!
+!--            Determine vertical index of topography top at grid point (j,i)
+               k_wall = MAXLOC(                                                &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr2 )   &
+                          ), DIM = 1                                           &
+                              ) - 1
+               DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
                ENDDO
@@ -4510,5 +5317,5 @@
 
 
-    SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
+    SUBROUTINE pmci_extrap_ifoutflow_lr( f, edge, var )
 !
 !--    After the interpolation of ghost-node values for the child-domain
@@ -4524,12 +5331,12 @@
        CHARACTER(LEN=1), INTENT(IN) ::  var    !:
 
-       INTEGER(iwp) ::  i     !:
-       INTEGER(iwp) ::  ib    !:
-       INTEGER(iwp) ::  ibgp  !:
-       INTEGER(iwp) ::  ied   !:
-       INTEGER(iwp) ::  j     !:
-       INTEGER(iwp) ::  k     !:
-      
-       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
+       INTEGER(iwp) ::  flag_nr !: Number of flag array to mask topography on respective u/v/w or s grid
+       INTEGER(iwp) ::  i       !:
+       INTEGER(iwp) ::  ib      !:
+       INTEGER(iwp) ::  ibgp    !:
+       INTEGER(iwp) ::  ied     !:
+       INTEGER(iwp) ::  j       !:
+       INTEGER(iwp) ::  k       !:
+       INTEGER(iwp) ::  k_wall  !:
 
        REAL(wp) ::  outnor    !:
@@ -4557,7 +5364,25 @@
           outnor = 1.0_wp
        ENDIF
+!
+!--    Determine number of flag array to be used to mask topography
+       IF ( var == 'u' )  THEN
+          flag_nr  = 14
+       ELSEIF ( var == 'v' )  THEN
+          flag_nr  = 16
+       ELSEIF ( var == 'w' )  THEN
+          flag_nr  = 18
+       ELSE
+          flag_nr  = 12
+       ENDIF
 
        DO  j = nys, nyn+1
-          DO  k = kb(j,i), nzt+1
+!
+!--       Determine vertical index of topography top at grid point (j,i)
+          k_wall = MAXLOC(                                                     &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr )    &
+                          ), DIM = 1                                           &
+                         ) - 1
+          DO  k = k_wall, nzt+1
              vdotnor = outnor * u(k,j,ied)
 !
@@ -4568,5 +5393,5 @@
           ENDDO
           IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
-             f(kb(j,i),j,i) = 0.0_wp
+             f(k_wall,j,i) = 0.0_wp
           ENDIF
        ENDDO
@@ -4588,5 +5413,5 @@
 
 
-    SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
+    SUBROUTINE pmci_extrap_ifoutflow_sn( f, edge, var )
 !
 !--    After  the interpolation of ghost-node values for the child-domain
@@ -4601,4 +5426,5 @@
        CHARACTER(LEN=1), INTENT(IN) ::  var    !:
       
+       INTEGER(iwp) ::  flag_nr   !: Number of flag array to mask topography on respective u/v/w or s grid
        INTEGER(iwp) ::  i         !:
        INTEGER(iwp) ::  j         !:
@@ -4607,6 +5433,5 @@
        INTEGER(iwp) ::  jed       !:
        INTEGER(iwp) ::  k         !:
-
-       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
+       INTEGER(iwp) ::  k_wall    !:
 
        REAL(wp)     ::  outnor    !:
@@ -4635,6 +5460,25 @@
        ENDIF
 
+!
+!--    Determine number of flag array to be used to mask topography
+       IF ( var == 'u' )  THEN
+          flag_nr  = 14
+       ELSEIF ( var == 'v' )  THEN
+          flag_nr  = 16
+       ELSEIF ( var == 'w' )  THEN
+          flag_nr  = 18
+       ELSE
+          flag_nr  = 12
+       ENDIF
+
        DO  i = nxl, nxr+1
-          DO  k = kb(j,i), nzt+1
+!
+!--       Determine vertical index of topography top at grid point (j,i)
+          k_wall = MAXLOC(                                                     &
+                     MERGE( 1, 0,                                              &
+                            BTEST( wall_flags_0(nzb:nzb_max,j,i), flag_nr )    &
+                          ), DIM = 1                                           &
+                         ) - 1
+          DO  k = k_wall, nzt+1
              vdotnor = outnor * v(k,jed,i)
 !
@@ -4645,5 +5489,5 @@
           ENDDO
           IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
-             f(kb(j,i),j,i) = 0.0_wp
+             f(k_wall,j,i) = 0.0_wp
           ENDIF
        ENDDO
