Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 2316)
+++ /palm/trunk/SOURCE/Makefile	(revision 2317)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# Add further dependencies on surface_mod
 # 
 # Former revisions:
@@ -521,5 +521,5 @@
 	lpm_splitting.o mod_kinds.o mod_particle_attributes.o
 lpm_advec.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o surface_mod.o
-lpm_boundary_conds.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o
+lpm_boundary_conds.o: modules.o cpulog_mod.o mod_kinds.o mod_particle_attributes.o surface_mod.o
 lpm_calc_liquid_water_content.o: modules.o cpulog_mod.o mod_kinds.o \
 	mod_particle_attributes.o
@@ -590,5 +590,5 @@
 	subsidence_mod.o surface_mod.o user_actions.o wind_turbine_model_mod.o
 progress_bar_mod.o: modules.o mod_kinds.o
-radiation_model_mod.o : modules.o mod_particle_attributes.o microphysics_mod.o
+radiation_model_mod.o: modules.o mod_particle_attributes.o microphysics_mod.o surface_mod.o
 random_function_mod.o: mod_kinds.o
 random_gauss.o: mod_kinds.o random_function_mod.o random_generator_parallel_mod.o
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 2316)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -278,5 +278,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  init_bc
+        ONLY:  get_topography_top_index, init_bc
 
     IMPLICIT NONE
@@ -1713,25 +1713,23 @@
 !--    nxl to nxr and nys to nyn on south and right model boundary, hence, 
 !--    use intrinsic lbound and ubound functions to infer array bounds.
-       DO  i = lbound(zu_s_inner, 1), ubound(zu_s_inner, 1)
-          DO  j = lbound(zu_s_inner, 2), ubound(zu_s_inner, 2)
+       DO  i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1)
+          DO  j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2)
 !
 !--          Topography height on scalar grid. Therefore, determine index of 
-!--          upward-facing surface element on scalar grid (bit 12). 
-             zu_s_inner(i,j) = zu( MAXLOC( MERGE(                              &
-                                         1, 0, BTEST( wall_flags_0(:,j,i), 12 )&
-                                                ), DIM = 1                     &
-                                         ) - 1                                 &
-                                 )
+!--          upward-facing surface element on scalar grid. 
+             zu_s_inner(i,j) = zu( get_topography_top_index( j, i, 's' ) )
+
+             write(9,*)  get_topography_top_index( j, i, 's' ), MAXLOC(                                                    &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(:,j,i), 12 )  &
+                               ), DIM = 1                                      &
+                          ) - 1
 !
 !--          Topography height on w grid. Therefore, determine index of 
-!--          upward-facing surface element on w grid (bit 18).
-             zw_w_inner(i,j) = zw( MAXLOC( MERGE(                              &
-                                         1, 0, BTEST( wall_flags_0(:,j,i), 18 )&
-                                                ), DIM = 1                     &
-                                         ) - 1                                 &
-                                 )
+!--          upward-facing surface element on w grid.
+             zw_w_inner(i,j) = zw( get_topography_top_index( j, i, 's' ) )
           ENDDO
        ENDDO
-
+       flush(9)
 
     ENDIF
Index: /palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_advec.f90	(revision 2316)
+++ /palm/trunk/SOURCE/lpm_advec.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -126,5 +126,5 @@
         
     USE indices,                                                               &
-        ONLY:  nzb, nzb_max, nzt, wall_flags_0
+        ONLY:  nzb, nzt
         
     USE kinds
@@ -140,5 +140,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -284,9 +284,5 @@
 !
 !--       Determine vertical index of topography top
-          k_wall = MAXLOC(                                                     &
-                       MERGE( 1, 0,                                            &
-                              BTEST( wall_flags_0(nzb:nzb_max,jlog,ilog), 12 ) &
-                            ), DIM = 1                                         &
-                         ) - 1
+          k_wall = get_topography_top_index( jlog,ilog, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -385,9 +381,5 @@
 !
 !--       Determine vertical index of topography top
-          k_wall = MAXLOC(                                                     &
-                       MERGE( 1, 0,                                            &
-                              BTEST( wall_flags_0(nzb:nzb_max,jlog,ilog), 12 ) &
-                            ), DIM = 1                                         &
-                         ) - 1
+          k_wall = get_topography_top_index( jlog,ilog, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -677,9 +669,5 @@
 !
 !--          Determine vertical index of topography top at (j,i)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
 !
 !--          To do: Reconsider order of computations in order to avoid 
@@ -700,9 +688,6 @@
 !
 !--          Determine vertical index of topography top at (j+1,i)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j+1,i), 12 )  &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j+1, i, 's' )
+
              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -720,9 +705,6 @@
 !
 !--          Determine vertical index of topography top at (j,i)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -740,9 +722,5 @@
 !
 !--          Determine vertical index of topography top at (j+1,i)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j+1,i), 12 )  &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j+1, i, 's' )
              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -760,9 +738,5 @@
 !
 !--          Determine vertical index of topography top at (j,i+1)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i+1), 12 )  &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i+1, 's' )
              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -780,9 +754,6 @@
 !
 !--          Determine vertical index of topography top at (j+1,i+1)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j+1,i+1), 12 )&
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j+1, i+1, 's' )
+
              IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -800,9 +771,6 @@
 !
 !--          Determine vertical index of topography top at (j,i+1)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i+1), 12 )  &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i+1, 's' )
+
              IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
                 num_gp = num_gp + 1
@@ -820,9 +788,6 @@
 !
 !--          Determine vertical index of topography top at (j+1,i+1)
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j+1,i+1), 12 )&
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j+1, i+1, 's' )
+
              IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
                 num_gp = num_gp + 1
Index: /palm/trunk/SOURCE/lpm_boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 2316)
+++ /palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -98,5 +98,5 @@
 
     USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys, nz, nzb, nzb_max, wall_flags_0
+        ONLY:  nxl, nxr, nyn, nys, nz, nzb
 
     USE kinds
@@ -108,4 +108,7 @@
 
     USE pegrid
+
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index
 
     IMPLICIT NONE
@@ -289,24 +292,8 @@
 !--       The construct of MERGE and BTEST is used to determine the topography-
 !--       top index (former nzb_s_inner). 
-          zwall1 = zw( MAXLOC(                                                 &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j2,i2), 12 )  &
-                               ), DIM = 1                                      &
-                             ) - 1 )                                             
-          zwall2 = zw( MAXLOC(                                                 &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j1,i1), 12 )  &
-                               ), DIM = 1                                      &
-                             ) - 1 )  
-          zwall3 = zw( MAXLOC(                                                 &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j1,i2), 12 )  &
-                               ), DIM = 1                                      &
-                             ) - 1 )  
-          zwall4 = zw( MAXLOC(                                                 &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j2,i1), 12 )  &
-                               ), DIM = 1                                      &
-                             ) - 1 )  
+          zwall1 = zw( get_topography_top_index( j2, i2, 's' ) )                                             
+          zwall2 = zw( get_topography_top_index( j1, i1, 's' ) ) 
+          zwall3 = zw( get_topography_top_index( j1, i2, 's' ) )
+          zwall4 = zw( get_topography_top_index( j2, i1, 's' ) )
 !
 !--       Initialize flags to check if particle reflection is necessary
@@ -493,9 +480,5 @@
 !--             necessarily exactly match the wall location due to rounding 
 !--             errors. At first, determine index of topography top at (j3,i3)  
-                k_wall = MAXLOC(                                               &
-                            MERGE( 1, 0,                                       &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )  &
-                                 ), DIM = 1                                    &
-                               ) - 1  
+                k_wall = get_topography_top_index( j3, i3, 's' ) 
                 IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
                      pos_z <= zw(k_wall)             .AND.                     &
@@ -534,9 +517,5 @@
 !--             necessary, carry out reflection. At first, determine index of
 !--             topography top at (j3,i3)  
-                k_wall = MAXLOC(                                               &
-                            MERGE( 1, 0,                                       &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )  &
-                                 ), DIM = 1                                    &
-                               ) - 1  
+                k_wall = get_topography_top_index( j3, i3, 's' ) 
                 IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
                      pos_z <= zw(k_wall)             .AND.                     &
@@ -565,9 +544,5 @@
 !--                Determine index of topography top at (j3,i3) and chick if
 !--                particle is below.  
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                   BTEST( wall_flags_0(nzb:nzb_max,j3,i3), 12 )&
-                               ), DIM = 1                                      &
-                                  ) - 1  
+                   k_wall = get_topography_top_index( j3, i3, 's' ) 
                    IF ( pos_z - zw(k_wall) < eps ) THEN
   
Index: /palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- /palm/trunk/SOURCE/lpm_init.f90	(revision 2316)
+++ /palm/trunk/SOURCE/lpm_init.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Get topography top index via Function call 
 !
 ! Former revisions:
@@ -178,5 +178,5 @@
     USE indices,                                                               &
         ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb,    &
-               nzb_max, nzt, wall_flags_0
+               nzt, wall_flags_0
 
     USE kinds
@@ -214,5 +214,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -725,11 +725,6 @@
 !
 !--                            Determine surface level. Therefore, check for
-!--                            upward-facing wall on w-grid. MAXLOC will return
-!--                            the index of the lowest upward-facing wall.
-                               k_surf = MAXLOC(                                &
-                                             MERGE( 1, 0,                      &
-                                   BTEST( wall_flags_0(nzb:nzb_max,jp,ip), 18 )&
-                                                  ), DIM = 1                   &
-                                              ) - 1
+!--                            upward-facing wall on w-grid. 
+                               k_surf = get_topography_top_index( j, i, 'w' )
 
                                IF ( seed_follows_topography )  THEN
Index: /palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- /palm/trunk/SOURCE/microphysics_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/microphysics_mod.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Get topography top index via Function call 
 !
 ! Former revisions:
@@ -445,5 +445,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -787,5 +787,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -923,7 +923,10 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
+
+       USE surface_mod,                                                        &
+           ONLY:  get_topography_top_index
 
 
@@ -942,9 +945,5 @@
 !
 !--          Determine vertical index of topography top
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
              DO  k = nzb+1, nzt
 !
@@ -994,5 +993,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1093,5 +1092,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1170,5 +1169,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1331,5 +1330,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1443,5 +1442,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -2329,7 +2328,10 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_max, nzt, wall_flags_0
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
+
+       USE surface_mod,                                                        &
+           ONLY:  get_topography_top_index
 
 
@@ -2346,9 +2348,5 @@
 !
 !--    Determine vertical index of topography top
-       k_wall = MAXLOC(                                                        &
-                        MERGE( 1, 0,                                           &
-                               BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )      &
-                             ), DIM = 1                                        &
-                      ) - 1
+       k_wall = get_topography_top_index( j, i, 's' )
        DO  k = nzb+1, nzt
 !
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -128,7 +128,10 @@
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
-               nz, nzb, nzb_max, nzt, wall_flags_0
+               nz, nzb, nzt
 
     USE kinds
+
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index
 
 
@@ -1016,9 +1019,5 @@
 !
 !--                Determine topography-top index on u-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 'u' )
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1084,9 +1083,6 @@
 !
 !--                Determine topography-top index on v-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 'v' )
+
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1152,9 +1148,6 @@
 !
 !--                Determine topography-top index on w-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 'w' )
+
                    DO  k = k_wall+1, k_wall+pch_index-1
 
@@ -1207,9 +1200,6 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 's' )
+
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1227,9 +1217,6 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 's' )
+
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1260,9 +1247,6 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 's' )
+
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1292,9 +1276,6 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = MAXLOC(                                            &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                                  ) - 1
+                   k_wall = get_topography_top_index( j, i, 's' )
+
                    DO  k = k_wall+1, k_wall+pch_index
 
@@ -1390,9 +1371,6 @@
 !
 !--          Determine topography-top index on u-grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 'u' )
+
              DO  k = k_wall+1, k_wall+pch_index
 
@@ -1453,9 +1431,6 @@
 !
 !--          Determine topography-top index on v-grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 'v' )
+
              DO  k = k_wall+1, k_wall+pch_index
 
@@ -1516,9 +1491,6 @@
 !
 !--          Determine topography-top index on w-grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 'w' )
+
              DO  k = k_wall+1, k_wall+pch_index-1
 
@@ -1566,9 +1538,6 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              DO  k = k_wall+1, k_wall+pch_index
                 kk = k - k_wall  !- lad arrays are defined flat
@@ -1582,9 +1551,6 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              DO  k = k_wall+1, k_wall+pch_index
 
@@ -1611,9 +1577,6 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              DO  k = k_wall+1, k_wall+pch_index
 
@@ -1640,9 +1603,6 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              DO  k = k_wall+1, k_wall+pch_index
 
Index: /palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2317)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -192,5 +192,5 @@
     USE indices,                                                                &
         ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg,  &
-               nysv, nz, nzb, nzb_max, nzt, wall_flags_0
+               nysv, nz, nzb, nzt, wall_flags_0
 
     USE kinds
@@ -233,5 +233,5 @@
 
     USE surface_mod,                                                            &
-        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -1394,33 +1394,17 @@
 !--             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          )
+                                         get_topography_top_index( j, i, 's' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'u' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'v' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'w' ) )
              ENDDO
           ENDDO
@@ -1436,33 +1420,17 @@
 !--             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          )
+                                         get_topography_top_index( j, i, 's' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'u' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'v' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'w' ) )
           ENDDO
           nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
@@ -1477,33 +1445,17 @@
 !--             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          )
+                                         get_topography_top_index( j, i, 's' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'u' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'v' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'w' ) )
              ENDDO
           ENDDO
@@ -1519,33 +1471,17 @@
 !--             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          )
+                                         get_topography_top_index( j, i, 's' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'u' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'v' ) )
 !
 !--             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          )
+                                         get_topography_top_index( j, i, 'w' ) )
           ENDDO
           nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
@@ -1603,8 +1539,5 @@
 !--          is part of the surfacetypes now. Set default roughness instead. 
 !--          Determine topography top index on u-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1622,8 +1555,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1669,8 +1599,5 @@
 !--          to the present surface tpye. 
 !--          Determine topography top index on u-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1688,8 +1615,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1732,8 +1656,5 @@
 !
 !--          Determine topography top index on u-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1751,8 +1672,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1794,9 +1712,6 @@
              j   = nyn + 1
 !
-!--          Determine topography top index on v-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+!--          Determine topography top index on u-grid
+             kb  = get_topography_top_index( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1814,8 +1729,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = MAXLOC( MERGE( 1, 0,                                        &
-                          BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) &
-                                ), DIM = 1                                     &
-                         ) - 1
+             kb  = get_topography_top_index( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1856,35 +1768,13 @@
 
              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
+!
+!--             Determine lowest grid on outer grids for u and w. 
+                k_wall_u_ji   = get_topography_top_index( j,   0, 'u_out' )
+                k_wall_u_ji_p = get_topography_top_index( j+1, 0, 'u_out' )
+                k_wall_u_ji_m = get_topography_top_index( j-1, 0, 'u_out' )
+
+                k_wall_w_ji   = get_topography_top_index( j,   -1, 'w_out' )
+                k_wall_w_ji_p = get_topography_top_index( j+1, -1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index( j-1, -1, 'w_out' )
 
                 DO  k = nzb, nzt_topo_nestbc_l
@@ -1971,36 +1861,14 @@
 
              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
+!
+!--             Determine lowest grid on outer grids for u and w. 
+                k_wall_u_ji   = get_topography_top_index( j,   i, 'u_out' )
+                k_wall_u_ji_p = get_topography_top_index( j+1, i, 'u_out' )
+                k_wall_u_ji_m = get_topography_top_index( j-1, i, 'u_out' )
+
+                k_wall_w_ji   = get_topography_top_index( j,   i, 'w_out' )
+                k_wall_w_ji_p = get_topography_top_index( j+1, i, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index( j-1, i, 'w_out' )
+
                 DO  k = nzb, nzt_topo_nestbc_r
 !
@@ -2081,36 +1949,14 @@
 
              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
+!
+!--             Determine lowest grid on outer grids for v and w. 
+                k_wall_v_ji   = get_topography_top_index( 0, i,   'v_out' )
+                k_wall_v_ji_p = get_topography_top_index( 0, i+1, 'v_out' )
+                k_wall_v_ji_m = get_topography_top_index( 0, i-1, 'v_out' )
+
+                k_wall_w_ji   = get_topography_top_index( -1, i,   'w_out' )
+                k_wall_w_ji_p = get_topography_top_index( -1, i+1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index( -1, i-1, 'w_out' )
+
                 DO  k = nzb, nzt_topo_nestbc_s
 !
@@ -2195,36 +2041,14 @@
 
              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
+!
+!--             Determine lowest grid on outer grids for v and w. 
+                k_wall_v_ji   = get_topography_top_index( j, i,   'v_out' )
+                k_wall_v_ji_p = get_topography_top_index( j, i+1, 'v_out' )
+                k_wall_v_ji_m = get_topography_top_index( j, i-1, 'v_out' )
+
+                k_wall_w_ji   = get_topography_top_index( j, i,   'w_out' )
+                k_wall_w_ji_p = get_topography_top_index( j, i+1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index( j, i-1, 'w_out' )
+
                 DO  k = nzb, nzt_topo_nestbc_n
 !
@@ -2831,9 +2655,5 @@
           i = nxl - 1
           DO  j = nysg, nyng
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2856,9 +2676,5 @@
           i = nxr + 1
           DO  j = nysg, nyng
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2881,9 +2697,5 @@
           j = nys - 1
           DO  i = nxlg, nxrg
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2906,9 +2718,6 @@
           j = nyn + 1
           DO  i = nxlg, nxrg
-             k_wall = MAXLOC(                                                  &
-                          MERGE( 1, 0,                                         &
-                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
-                               ), DIM = 1                                      &
-                            ) - 1
+             k_wall = get_topography_top_index( j, i, 's' )
+
              DO  k = k_wall + 1, nzt
 
@@ -2932,9 +2741,5 @@
 !
 !--          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
+             k_wall = get_topography_top_index( j, i, 's' )
 
              kc     = kco(k+1)
@@ -3315,6 +3120,4 @@
        INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
 
-       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       !:
@@ -3380,19 +3183,5 @@
           ENDIF
        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.
@@ -3426,9 +3215,5 @@
 !
 !--             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_wall = get_topography_top_index( j, i, TRIM ( var ) )
 !
 !--             kbc is the first coarse-grid point above the surface
@@ -3461,9 +3246,5 @@
 !
 !--              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_wall = get_topography_top_index( j, i, 'w' )
 
                  f(k_wall,j,i) = 0.0_wp
@@ -4563,6 +4344,4 @@
       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) ::  ib       !:
@@ -4611,19 +4390,4 @@
          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
@@ -4653,9 +4417,5 @@
 !
 !--         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_wall = get_topography_top_index( j, i, TRIM ( var ) )
 
             k = k_wall+1
@@ -4683,9 +4443,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, TRIM ( var ) )
+
                DO  k = k_wall+1, nzt_topo_nestbc
                   IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) == 0 ) )  THEN
@@ -4711,9 +4468,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, TRIM ( var ) )
+
                k = k_wall + 1
                IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) /= 0 ) )  THEN
@@ -4746,9 +4500,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, 's' )
+
                DO  k = k_wall, nzt + 1
                   f(k,j,i) = tkefactor_l(k,j) * f(k,j,i)
@@ -4759,9 +4510,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, 's' )
+
                DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
@@ -4821,6 +4569,4 @@
       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    !:
@@ -4867,28 +4613,10 @@
       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
 !
 !--      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_wall = get_topography_top_index( j, i, TRIM( var ) )
+
          DO  k = k_wall, nzt+1
             l = ic(i)
@@ -4916,9 +4644,5 @@
 !
 !--         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_wall = get_topography_top_index( j, i, TRIM( var ) )
 
             k = k_wall + 1
@@ -4945,9 +4669,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, TRIM( var ) )
+
                DO  k = k_wall, nzt_topo_nestbc
 !
@@ -4975,9 +4696,6 @@
 !
 !--            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_wall = get_topography_top_index( j, i, TRIM( var ) )
+
                k = k_wall + 1
                IF ( ( logc(2,k,i) /= 0 )  .AND.  ( logc(1,k,i) /= 0 ) )  THEN
@@ -5010,9 +4728,5 @@
 !
 !--            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_wall = get_topography_top_index( j, i, 's' )
                DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
@@ -5023,9 +4737,5 @@
 !
 !--            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_wall = get_topography_top_index( j, i, 's' )
                DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
@@ -5153,5 +4863,4 @@
        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) ::  ib      !:
@@ -5186,24 +4895,11 @@
           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
 !
 !--       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
+          k_wall = get_topography_top_index( j, i, TRIM( var ) )
+
           DO  k = k_wall, nzt+1
              vdotnor = outnor * u(k,j,ied)
@@ -5248,5 +4944,4 @@
        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         !:
@@ -5282,24 +4977,10 @@
        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
 !
 !--       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
+          k_wall = get_topography_top_index( j, i, TRIM( var ) )
+
           DO  k = k_wall, nzt+1
              vdotnor = outnor * v(k,jed,i)
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -169,6 +169,5 @@
 
     USE indices,                                                               &
-        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
-               wall_flags_0
+        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
 
     USE kinds
@@ -200,6 +199,6 @@
         ONLY:  rrtmg_sw
 #endif
-
-
+    USE surface_mod,                                                           &
+        ONLY:  get_topography_top_index
 
     IMPLICIT NONE
@@ -1186,9 +1185,5 @@
 !
 !--          Obtain vertical index of topography top 
-             k = MAXLOC(                                                       &
-                        MERGE( 1, 0,                                           &
-                               BTEST( wall_flags_0(:,j,i), 12 )                &
-                             ), DIM = 1                                        &
-                       ) - 1
+             k = get_topography_top_index( j, i, 's' )
 
              exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
@@ -1243,9 +1238,5 @@
 !--          Obtain vertical index of topography top. So far it is identical to
 !--          nzb.
-             k = MAXLOC(                                                       &
-                        MERGE( 1, 0,                                           &
-                               BTEST( wall_flags_0(:,j,i), 12 )                &
-                             ), DIM = 1                                        &
-                       ) - 1
+             k = get_topography_top_index( j, i, 's' )
 
              rad_net(j,i)      = net_radiation
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 2317)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! New function to obtain topography top index. 
 ! 
 ! Former revisions:
@@ -269,4 +269,33 @@
 
     PRIVATE
+
+    INTERFACE get_topography_top_index
+       MODULE PROCEDURE get_topography_top_index
+    END  INTERFACE get_topography_top_index
+
+    INTERFACE init_bc
+       MODULE PROCEDURE init_bc
+    END INTERFACE init_bc
+
+    INTERFACE init_surfaces
+       MODULE PROCEDURE init_surfaces
+    END INTERFACE init_surfaces
+
+    INTERFACE init_surface_arrays
+       MODULE PROCEDURE init_surface_arrays
+    END INTERFACE init_surface_arrays
+
+    INTERFACE surface_read_restart_data
+       MODULE PROCEDURE surface_read_restart_data
+    END INTERFACE surface_read_restart_data
+
+    INTERFACE surface_write_restart_data
+       MODULE PROCEDURE surface_write_restart_data
+    END INTERFACE surface_write_restart_data
+
+    INTERFACE surface_last_actions
+       MODULE PROCEDURE surface_last_actions
+    END INTERFACE surface_last_actions
+
 !
 !-- Public variables
@@ -274,32 +303,8 @@
            surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type
 !
-!-- Public subroutines
-    PUBLIC init_bc, init_surfaces, init_surface_arrays,                        &
-           surface_read_restart_data, surface_write_restart_data,              &
-           surface_last_actions
-
-    INTERFACE init_bc
-       MODULE PROCEDURE init_bc
-    END INTERFACE init_bc
-
-    INTERFACE init_surfaces
-       MODULE PROCEDURE init_surfaces
-    END INTERFACE init_surfaces
-
-    INTERFACE init_surface_arrays
-       MODULE PROCEDURE init_surface_arrays
-    END INTERFACE init_surface_arrays
-
-    INTERFACE surface_read_restart_data
-       MODULE PROCEDURE surface_read_restart_data
-    END INTERFACE surface_read_restart_data
-
-    INTERFACE surface_write_restart_data
-       MODULE PROCEDURE surface_write_restart_data
-    END INTERFACE surface_write_restart_data
-
-    INTERFACE surface_last_actions
-       MODULE PROCEDURE surface_last_actions
-    END INTERFACE surface_last_actions
+!-- Public subroutines and functions
+    PUBLIC get_topography_top_index, init_bc, init_surfaces,                   &
+           init_surface_arrays, surface_read_restart_data,                     &
+           surface_write_restart_data, surface_last_actions
 
 
@@ -1642,4 +1647,57 @@
 
     END SUBROUTINE init_surfaces
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Determines topography-top index at given (j,i)-position.  
+!------------------------------------------------------------------------------!
+    FUNCTION get_topography_top_index( j, i, grid )
+
+       USE kinds
+
+       IMPLICIT NONE 
+
+       CHARACTER(LEN=*) ::  grid                      !< flag to distinquish between staggered grids
+       INTEGER(iwp)     ::  i                         !< grid index in x-dimension
+       INTEGER(iwp)     ::  ibit                      !< bit position where topography information is stored on respective grid
+       INTEGER(iwp)     ::  j                         !< grid index in y-dimension
+       INTEGER(iwp)     ::  get_topography_top_index  !< topography top index
+
+       SELECT CASE ( TRIM( grid ) )
+
+          CASE ( 's'     )
+             ibit = 12
+          CASE ( 'u'     )
+             ibit = 14
+          CASE ( 'v'     )
+             ibit = 16
+          CASE ( 'w'     )
+             ibit = 18
+          CASE ( 's_out' )
+             ibit = 24
+          CASE ( 'u_out' )
+             ibit = 26
+          CASE ( 'v_out' )
+             ibit = 27
+          CASE ( 'w_out' )
+             ibit = 28
+          CASE DEFAULT
+!
+!--          Set default to scalar grid
+             ibit = 12 
+
+       END SELECT
+
+       get_topography_top_index = MAXLOC(                                      &
+                                     MERGE( 1, 0,                              &
+                                            BTEST( wall_flags_0(:,j,i), ibit ) &
+                                          ), DIM = 1                           &
+                                        ) - 1
+
+       RETURN
+
+    END FUNCTION get_topography_top_index
 
 !------------------------------------------------------------------------------!
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2316)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2317)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Get topography top index via Function call 
 ! 
 ! Former revisions:
@@ -603,8 +603,6 @@
 !
 !--                 Find topography top index
-                    k_topo = MAXLOC( MERGE( 1, 0,                              &
-                                             BTEST( wall_flags_0(:,j,i), 12 )  &
-                                           ), DIM = 1                          &
-                                   ) - 1
+                    k_topo = get_topography_top_index( j, i, 's' )
+
                     DO k = nzt+1, 0, -1
                         IF ( lad_s(k,j,i) /= 0.0_wp )  THEN
@@ -706,12 +704,6 @@
                  DO  j = ijdb(3,ids), ijdb(4,ids)
 
-                    k_topo =   MAXLOC( MERGE( 1, 0,                            &
-                                             BTEST( wall_flags_0(:,j,i), 12 )  &
-                                            ), DIM = 1                         &
-                                     ) - 1
-                    k_topo2 =  MAXLOC( MERGE( 1, 0,                            &
-                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
-                                            ), DIM = 1                         &
-                                     ) - 1 
+                    k_topo  = get_topography_top_index( j, i, 's' )
+                    k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
 
                     k = nzut - MAX( k_topo, k_topo2 )
@@ -732,8 +724,6 @@
 !
 !--                 Find topography top index
-                    k_topo = MAXLOC( MERGE( 1, 0,                              &
-                                             BTEST( wall_flags_0(:,j,i), 12 )  &
-                                          ), DIM = 1                           &
-                                   ) - 1
+                    k_topo = get_topography_top_index( j, i, 's' )
+
                     DO k = k_topo + 1, pct(j,i)
                         ipcgb = ipcgb + 1
@@ -813,12 +803,7 @@
                 DO i = ijdb(1,ids), ijdb(2,ids)
                     DO j = ijdb(3,ids), ijdb(4,ids)
-                        k_topo =   MAXLOC( MERGE( 1, 0,                        &
-                                             BTEST( wall_flags_0(:,j,i), 12 )  &
-                                                ), DIM = 1                     &
-                                         ) - 1
-                        k_topo2 =  MAXLOC( MERGE( 1, 0,                        &
-                                             BTEST( wall_flags_0(:,j-jdir(ids),i-idir(ids)), 12 )  &
-                                                ), DIM = 1                     &
-                                         ) - 1 
+                        k_topo  = get_topography_top_index( j, i, 's' )
+                        k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
+
                         DO k = MAX(k_topo,k_topo2)+1, nzut
                             isurf = isurf + 1
@@ -1929,9 +1914,5 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
-                    k = MAXLOC(                                                &
-                                MERGE( 1, 0,                                   &
-                                       BTEST( wall_flags_0(:,j,i), 12 )        &
-                                     ), DIM = 1                                &
-                              ) - 1
+                    k = get_topography_top_index( j, i, 's' )
 
                     usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
@@ -3806,10 +3787,5 @@
 !
 !--             Following expression equals former kk = k - nzb_s_inner(j,i) 
-                kk = k - ( MAXLOC(                                             &
-                                MERGE( 1, 0,                                   &
-                                       BTEST( wall_flags_0(:,j,i), 12 )        &
-                                     ), DIM = 1                                &
-                               ) - 1                                           &
-                         )  !- lad arrays are defined flat
+                kk = k - ( get_topography_top_index( j, i, 's' ) )  !- lad arrays are defined flat
                 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
                     * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
