Index: palm/trunk/SOURCE/date_and_time_mod.f90
===================================================================
--- palm/trunk/SOURCE/date_and_time_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/date_and_time_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in definition of d_seconds_year.
 ! 
 ! Former revisions:
@@ -56,7 +56,7 @@
     REAL(wp) ::  time_utc_init = 43200.0_wp   !< UTC time at model start
 
-    REAL(wp), PARAMETER ::  d_hours_day    = 0.0416666666667_wp    !< inverse of hours per day (1/24)
-    REAL(wp), PARAMETER ::  d_seconds_hour = 0.000277777777778_wp   !< inverse of seconds per hour (1/3600)
-    REAL(wp), PARAMETER ::  d_seconds_year = 31536000.0_wp          !< inverse of the seconds per year (1/(365*86400))
+    REAL(wp), PARAMETER ::  d_hours_day    = 1.0_wp / 24.0_wp       !< inverse of hours per day (1/24)
+    REAL(wp), PARAMETER ::  d_seconds_hour = 1.0_wp / 3600.0_wp     !< inverse of seconds per hour (1/3600)
+    REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
     
     SAVE
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 2697)
+++ palm/trunk/SOURCE/header.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -410,5 +410,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, get_topography_top_index
+        ONLY:  surf_def_h, get_topography_top_index_ji
 
     USE synthetic_turbulence_generator_mod,                                    &
@@ -993,5 +993,5 @@
        WRITE( io, 280 ) 
        IF ( turbulent_inflow )  THEN
-          WRITE( io, 281 )  zu( get_topography_top_index( 0, 0, 's' ) )
+          WRITE( io, 281 )  zu( get_topography_top_index_ji( 0, 0, 's' ) )
        ENDIF
        IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 2697)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -466,5 +466,5 @@
     USE surface_mod,                                                           &
         ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
-                surf_usm_h, get_topography_top_index
+                surf_usm_h, get_topography_top_index_ji
    
     USE transpose_indices
@@ -1580,8 +1580,8 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                nz_u_shift = get_topography_top_index( j, i, 'u' )
-                nz_v_shift = get_topography_top_index( j, i, 'v' )
-                nz_w_shift = get_topography_top_index( j, i, 'w' )
-                nz_s_shift = get_topography_top_index( j, i, 's' )
+                nz_u_shift = get_topography_top_index_ji( j, i, 'u' )
+                nz_v_shift = get_topography_top_index_ji( j, i, 'v' )
+                nz_w_shift = get_topography_top_index_ji( j, i, 'w' )
+                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
 
                 u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)                
@@ -1625,8 +1625,8 @@
           IF ( complex_terrain )  THEN
              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
-                nz_u_shift_l = get_topography_top_index( 0, 0, 'u' )
-                nz_v_shift_l = get_topography_top_index( 0, 0, 'v' )
-                nz_w_shift_l = get_topography_top_index( 0, 0, 'w' )
-                nz_s_shift_l = get_topography_top_index( 0, 0, 's' )
+                nz_u_shift_l = get_topography_top_index_ji( 0, 0, 'u' )
+                nz_v_shift_l = get_topography_top_index_ji( 0, 0, 'v' )
+                nz_w_shift_l = get_topography_top_index_ji( 0, 0, 'w' )
+                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
              ELSE
                 nz_u_shift_l = 0
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 2697)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -298,5 +298,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, init_bc
+        ONLY:  get_topography_top_index, get_topography_top_index_ji, init_bc
 
     USE vertical_nesting_mod,                                                  &
@@ -552,9 +552,9 @@
 !--          Topography height on scalar grid. Therefore, determine index of 
 !--          upward-facing surface element on scalar grid. 
-             zu_s_inner(i,j) = zu( get_topography_top_index( j, i, 's' ) )
+             zu_s_inner(i,j) = zu( get_topography_top_index_ji( j, i, 's' ) )
 !
 !--          Topography height on w grid. Therefore, determine index of 
 !--          upward-facing surface element on w grid.
-             zw_w_inner(i,j) = zw( get_topography_top_index( j, i, 's' ) )
+             zw_w_inner(i,j) = zw( get_topography_top_index_ji( j, i, 's' ) )
           ENDDO
        ENDDO
@@ -1471,5 +1471,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index
+        ONLY:  get_topography_top_index, get_topography_top_index_ji
 
     IMPLICIT NONE
Index: palm/trunk/SOURCE/lpm.f90
===================================================================
--- palm/trunk/SOURCE/lpm.f90	(revision 2697)
+++ palm/trunk/SOURCE/lpm.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele)
 ! 
 ! Former revisions:
@@ -343,5 +343,5 @@
 !--             optimization is still possible.)
                 IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
-                   CALL lpm_boundary_conds( 'walls' )
+                   CALL  lpm_boundary_conds( 'walls', i, j, k )
                 ENDIF
 !
@@ -353,5 +353,5 @@
 !--             the top or bottom boundary and delete those particles, which are
 !--             older than allowed
-                CALL lpm_boundary_conds( 'bottom/top' )
+                CALL lpm_boundary_conds( 'bottom/top', i, j, k )
 !
 !---            If not all particles of the actual grid cell have reached the 
Index: palm/trunk/SOURCE/lpm_advec.f90
===================================================================
--- palm/trunk/SOURCE/lpm_advec.f90	(revision 2697)
+++ palm/trunk/SOURCE/lpm_advec.f90	(revision 2698)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! Particle interpolations at walls in case of SGS velocities revised and not
+! required parts are removed. (responsible Philipp Thiele)
+!
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -149,5 +152,5 @@
         
     USE indices,                                                               &
-        ONLY:  nzb, nzt
+        ONLY:  nzb, nzt, wall_flags_0
         
     USE kinds
@@ -163,7 +166,9 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
+
+    LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
 
     INTEGER(iwp) ::  agp                         !< loop variable
@@ -302,5 +307,5 @@
 !--       above topography (Prandtl-layer height)
 !--       Determine vertical index of topography top
-          k_wall = get_topography_top_index( jp, ip, 's' )
+          k_wall = get_topography_top_index_ji( jp, ip, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -397,5 +402,5 @@
 !
 !--       Determine vertical index of topography top
-          k_wall = get_topography_top_index( jp,ip, 's' )
+          k_wall = get_topography_top_index_ji( jp,ip, 's' )
 
           IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -528,9 +533,31 @@
 !-- velocities
     IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
-
-       IF ( topography == 'flat' )  THEN
-
-          DO  nb = 0,7
-
+    
+       DO  nb = 0,7
+         
+          subbox_at_wall = .FALSE.         
+!
+!--       In case of topography check if subbox is adjacent to a wall
+          IF ( .NOT. topography == 'flat' ) THEN
+             i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
+             j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
+             k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) )
+             IF ( .NOT. BTEST(wall_flags_0(k,  jp, ip), 0) .OR.                &
+                  .NOT. BTEST(wall_flags_0(kp, j,  ip), 0) .OR.                &
+                  .NOT. BTEST(wall_flags_0(kp, jp, i ), 0) )                   &
+             THEN
+                subbox_at_wall = .TRUE.
+             ENDIF
+          ENDIF
+          IF ( subbox_at_wall ) THEN
+             e_int(start_index(nb):end_index(nb))     = e(kp,jp,ip) 
+             diss_int(start_index(nb):end_index(nb))  = diss(kp,jp,ip)
+             de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip)
+             de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip)
+             de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip)
+!
+!--          Set flag for stochastic equation.
+             term_1_2(start_index(nb):end_index(nb)) = 0.0_wp             
+          ELSE
              i = ip + block_offset(nb)%i_off
              j = jp + block_offset(nb)%j_off
@@ -600,10 +627,10 @@
                 ELSE
                    de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
-                                  ( gg - bb ) * de_dy(k+1,j,i+1) + &
-                                  ( gg - cc ) * de_dy(k+1,j+1,i) + &
-                                  ( gg - dd ) * de_dy(k+1,j+1,i+1) &
+                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
+                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
+                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
                                   ) / ( 3.0_wp * gg )
-                   de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
-                                              ( de_dy_int_u - de_dy_int_l )
+                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
+                                                 ( de_dy_int_u - de_dy_int_l )
                 ENDIF
 
@@ -638,5 +665,5 @@
                                ( gg - cc ) * diss(k,j+1,i) + &
                                ( gg - dd ) * diss(k,j+1,i+1) &
-                              ) / ( 3.0_wp * gg )
+                               ) / ( 3.0_wp * gg )
 
                 IF ( k == nzt )  THEN
@@ -649,5 +676,5 @@
                                  ) / ( 3.0_wp * gg )
                    diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
-                                           ( diss_int_u - diss_int_l )
+                                            ( diss_int_u - diss_int_l )
                 ENDIF
 
@@ -655,657 +682,7 @@
 !--             Set flag for stochastic equation.
                 term_1_2(n) = 1.0_wp
-
              ENDDO
-          ENDDO
-
-       ELSE ! non-flat topography, e.g., buildings 
-
-          DO  nb = 0, 7
-
-             DO  n = start_index(nb), end_index(nb)
-
-                i = particles(n)%x * ddx
-                j = particles(n)%y * ddy
-                k = ( zv(n) + dz * atmos_ocean_sign ) / dz  &
-                    + offset_ocean_nzt                      ! only exact if eq.dist
-!
-!--             In case that there are buildings it has to be determined 
-!--             how many of the gridpoints defining the particle box are 
-!--             situated within a building
-!--             gp_outside_of_building(1): i,j,k
-!--             gp_outside_of_building(2): i,j+1,k
-!--             gp_outside_of_building(3): i,j,k+1
-!--             gp_outside_of_building(4): i,j+1,k+1
-!--             gp_outside_of_building(5): i+1,j,k
-!--             gp_outside_of_building(6): i+1,j+1,k
-!--             gp_outside_of_building(7): i+1,j,k+1
-!--             gp_outside_of_building(8): i+1,j+1,k+1
-
-                gp_outside_of_building = 0
-                location = 0.0_wp
-                num_gp = 0
-
-!
-!--             Determine vertical index of topography top at (j,i)
-                k_wall = get_topography_top_index( j, i, 's' )
-!
-!--             To do: Reconsider order of computations in order to avoid 
-!--             unnecessary index calculations.
-                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN 
-                   num_gp = num_gp + 1
-                   gp_outside_of_building(1) = 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j,i)
-                   dissi(num_gp)  = diss(k,j,i)
-                   de_dxi(num_gp) = de_dx(k,j,i)
-                   de_dyi(num_gp) = de_dy(k,j,i)
-                   de_dzi(num_gp) = de_dz(k,j,i)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j+1,i)
-                k_wall = get_topography_top_index( j+1, i, 's' )
-
-                IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
-                   num_gp = num_gp + 1
-                   gp_outside_of_building(2) = 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j+1,i)
-                   dissi(num_gp)  = diss(k,j+1,i)
-                   de_dxi(num_gp) = de_dx(k,j+1,i)
-                   de_dyi(num_gp) = de_dy(k,j+1,i)
-                   de_dzi(num_gp) = de_dz(k,j+1,i)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j,i)
-                k_wall = get_topography_top_index( j, i, 's' )
-
-                IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
-                   num_gp = num_gp + 1
-                   gp_outside_of_building(3) = 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i)
-                   dissi(num_gp)  = diss(k+1,j,i)
-                   de_dxi(num_gp) = de_dx(k+1,j,i)
-                   de_dyi(num_gp) = de_dy(k+1,j,i)
-                   de_dzi(num_gp) = de_dz(k+1,j,i)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j+1,i)
-                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
-                   gp_outside_of_building(4) = 1
-                   location(num_gp,1) = i * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i)
-                   dissi(num_gp)  = diss(k+1,j+1,i)
-                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j,i+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
-                   gp_outside_of_building(5) = 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j,i+1)
-                   dissi(num_gp)  = diss(k,j,i+1)
-                   de_dxi(num_gp) = de_dx(k,j,i+1)
-                   de_dyi(num_gp) = de_dy(k,j,i+1)
-                   de_dzi(num_gp) = de_dz(k,j,i+1)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j+1,i+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
-                   gp_outside_of_building(6) = 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = k * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k,j+1,i+1)
-                   dissi(num_gp)  = diss(k,j+1,i+1)
-                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
-                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
-                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j,i+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
-                   gp_outside_of_building(7) = 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = j * dy
-                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j,i+1)
-                   dissi(num_gp)  = diss(k+1,j,i+1)
-                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
-                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
-                ENDIF
-
-!
-!--             Determine vertical index of topography top at (j+1,i+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
-                   gp_outside_of_building(8) = 1
-                   location(num_gp,1) = (i+1) * dx
-                   location(num_gp,2) = (j+1) * dy
-                   location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
-                   ei(num_gp)     = e(k+1,j+1,i+1)
-                   dissi(num_gp)  = diss(k+1,j+1,i+1)
-                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
-                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-                ENDIF
-!
-!--             If all gridpoints are situated outside of a building, then the
-!--             ordinary interpolation scheme can be used.
-                IF ( num_gp == 8 )  THEN
-
-                   x  = particles(n)%x - i * dx
-                   y  = particles(n)%y - j * dy
-                   aa = x**2          + y**2
-                   bb = ( dx - x )**2 + y**2
-                   cc = x**2          + ( dy - y )**2
-                   dd = ( dx - x )**2 + ( dy - y )**2
-                   gg = aa + bb + cc + dd
-     
-                   e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
-                             + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
-                             ) / ( 3.0_wp * gg )
-      
-                   IF ( k == nzt )  THEN
-                      e_int(n) = e_int_l
-                   ELSE
-                      e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
-                                  ( gg - bb ) * e(k+1,j,i+1) + &
-                                  ( gg - cc ) * e(k+1,j+1,i) + &
-                                  ( gg - dd ) * e(k+1,j+1,i+1) &
-                                ) / ( 3.0_wp * gg )
-                      e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
-                                          ( e_int_u - e_int_l )
-                   ENDIF
-! 
-!--                Needed to avoid NaN particle velocities (this might not be
-!--                required any more)
-                   IF ( e_int(n) <= 0.0_wp )  THEN
-                      e_int(n) = 1.0E-20_wp
-                   ENDIF
-!
-!--                Interpolate the TKE gradient along x (adopt incides i,j,k
-!--                and all position variables from above (TKE))
-                   de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
-                                   ( gg - bb ) * de_dx(k,j,i+1) + &
-                                   ( gg - cc ) * de_dx(k,j+1,i) + &
-                                   ( gg - dd ) * de_dx(k,j+1,i+1) &
-                                 ) / ( 3.0_wp * gg )
-
-                   IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
-                      de_dx_int(n) = de_dx_int_l
-                   ELSE
-                      de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
-                                      ( gg - bb ) * de_dx(k+1,j,i+1) + &
-                                      ( gg - cc ) * de_dx(k+1,j+1,i) + &
-                                      ( gg - dd ) * de_dx(k+1,j+1,i+1) &
-                                    ) / ( 3.0_wp * gg )
-                      de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
-                                              dz * ( de_dx_int_u - de_dx_int_l )
-                   ENDIF
-
-!
-!--                Interpolate the TKE gradient along y
-                   de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
-                                   ( gg - bb ) * de_dy(k,j,i+1) + &
-                                   ( gg - cc ) * de_dy(k,j+1,i) + &
-                                   ( gg - dd ) * de_dy(k,j+1,i+1) &
-                                 ) / ( 3.0_wp * gg )
-                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                      de_dy_int(n) = de_dy_int_l
-                   ELSE
-                      de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
-                                      ( gg - bb ) * de_dy(k+1,j,i+1) + &
-                                      ( gg - cc ) * de_dy(k+1,j+1,i) + &
-                                      ( gg - dd ) * de_dy(k+1,j+1,i+1) &
-                                    ) / ( 3.0_wp * gg )
-                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
-                                              dz * ( de_dy_int_u - de_dy_int_l )
-                   ENDIF
-
-!
-!--                Interpolate the TKE gradient along z
-                   IF ( zv(n) < 0.5_wp * dz )  THEN
-                      de_dz_int(n) = 0.0_wp
-                   ELSE
-                      de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
-                                      ( gg - bb ) * de_dz(k,j,i+1) + &
-                                      ( gg - cc ) * de_dz(k,j+1,i) + &
-                                      ( gg - dd ) * de_dz(k,j+1,i+1) &
-                                    ) / ( 3.0_wp * gg )
-    
-                      IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
-                         de_dz_int(n) = de_dz_int_l
-                      ELSE
-                         de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
-                                         ( gg - bb ) * de_dz(k+1,j,i+1) + &
-                                         ( gg - cc ) * de_dz(k+1,j+1,i) + &
-                                         ( gg - dd ) * de_dz(k+1,j+1,i+1) &
-                                       ) / ( 3.0_wp * gg )
-                         de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
-                                              dz * ( de_dz_int_u - de_dz_int_l )
-                      ENDIF
-                   ENDIF
-
-!
-!--                Interpolate the dissipation of TKE
-                   diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
-                                  ( gg - bb ) * diss(k,j,i+1) + &
-                                  ( gg - cc ) * diss(k,j+1,i) + &
-                                  ( gg - dd ) * diss(k,j+1,i+1) &
-                                ) / ( 3.0_wp * gg )
-    
-                   IF ( k == nzt )  THEN
-                      diss_int(n) = diss_int_l
-                   ELSE
-                      diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
-                                     ( gg - bb ) * diss(k+1,j,i+1) + &
-                                     ( gg - cc ) * diss(k+1,j+1,i) + &
-                                     ( gg - dd ) * diss(k+1,j+1,i+1) &
-                                   ) / ( 3.0_wp * gg )
-                      diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
-                                              ( diss_int_u - diss_int_l )
-                   ENDIF
-!
-!--                Set flag for stochastic equation.
-                   term_1_2(n) = 1.0_wp
-    
-                ELSE
-    
-!
-!--                If wall between gridpoint 1 and gridpoint 5, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(1) == 1  .AND. &
-                        gp_outside_of_building(5) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j,i)
-                      dissi(num_gp)  = diss(k,j,i)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k,j,i)
-                      de_dzi(num_gp) = de_dz(k,j,i)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(5) == 1  .AND. &
-                        gp_outside_of_building(1) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j,i+1)
-                      dissi(num_gp)  = diss(k,j,i+1)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k,j,i+1)
-                      de_dzi(num_gp) = de_dz(k,j,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 5 and gridpoint 6, then
-!--                then Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(5) == 1  .AND. &
-                        gp_outside_of_building(6) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j,i+1)
-                      dissi(num_gp)  = diss(k,j,i+1)
-                      de_dxi(num_gp) = de_dx(k,j,i+1)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k,j,i+1)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(6) == 1  .AND. &
-                        gp_outside_of_building(5) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j+1,i+1)
-                      dissi(num_gp)  = diss(k,j+1,i+1)
-                      de_dxi(num_gp) = de_dx(k,j+1,i+1)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 2 and gridpoint 6, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(2) == 1  .AND. &
-                        gp_outside_of_building(6) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j+1,i)
-                      dissi(num_gp)  = diss(k,j+1,i)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k,j+1,i)
-                      de_dzi(num_gp) = de_dz(k,j+1,i)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(6) == 1  .AND. &
-                        gp_outside_of_building(2) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j+1,i+1)
-                      dissi(num_gp)  = diss(k,j+1,i+1)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k,j+1,i+1)
-                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 1 and gridpoint 2, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(1) == 1  .AND. &
-                        gp_outside_of_building(2) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j,i)
-                      dissi(num_gp)  = diss(k,j,i)
-                      de_dxi(num_gp) = de_dx(k,j,i)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k,j,i)  
-                   ENDIF
-
-                   IF ( gp_outside_of_building(2) == 1  .AND. &
-                        gp_outside_of_building(1) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz - 0.5_wp * dz
-                      ei(num_gp)     = e(k,j+1,i)
-                      dissi(num_gp)  = diss(k,j+1,i)
-                      de_dxi(num_gp) = de_dx(k,j+1,i)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k,j+1,i)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 3 and gridpoint 7, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(3) == 1  .AND. &
-                        gp_outside_of_building(7) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j,i)
-                      dissi(num_gp)  = diss(k+1,j,i)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k+1,j,i)
-                      de_dzi(num_gp) = de_dz(k+1,j,i) 
-                   ENDIF
-
-                   IF ( gp_outside_of_building(7) == 1  .AND. &
-                        gp_outside_of_building(3) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j,i+1)
-                      dissi(num_gp)  = diss(k+1,j,i+1)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 7 and gridpoint 8, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(7) == 1  .AND. &
-                        gp_outside_of_building(8) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j,i+1)
-                      dissi(num_gp)  = diss(k+1,j,i+1)
-                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(8) == 1  .AND. &
-                        gp_outside_of_building(7) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j+1,i+1)
-                      dissi(num_gp)  = diss(k+1,j+1,i+1)
-                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 4 and gridpoint 8, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(4) == 1  .AND. &
-                        gp_outside_of_building(8) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j+1,i)
-                      dissi(num_gp)  = diss(k+1,j+1,i)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(8) == 1  .AND. &
-                        gp_outside_of_building(4) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx + 0.5_wp * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j+1,i+1)
-                      dissi(num_gp)  = diss(k+1,j+1,i+1)
-                      de_dxi(num_gp) = 0.0_wp
-                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 3 and gridpoint 4, then
-!--                Neumann boundary condition has to be applied
-                   IF ( gp_outside_of_building(3) == 1  .AND. &
-                        gp_outside_of_building(4) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j,i)
-                      dissi(num_gp)  = diss(k+1,j,i)
-                      de_dxi(num_gp) = de_dx(k+1,j,i)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k+1,j,i)
-                   ENDIF
-
-                   IF ( gp_outside_of_building(4) == 1  .AND. &
-                        gp_outside_of_building(3) == 0 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = j * dy + 0.5_wp * dy
-                      location(num_gp,3) = k * dz + 0.5_wp * dz
-                      ei(num_gp)     = e(k+1,j+1,i)
-                      dissi(num_gp)  = diss(k+1,j+1,i)
-                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                      de_dyi(num_gp) = 0.0_wp
-                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
-                   ENDIF
-
-!
-!--                If wall between gridpoint 1 and gridpoint 3, then
-!--                Neumann boundary condition has to be applied
-!--                (only one case as only building beneath is possible)
-                   IF ( gp_outside_of_building(1) == 0  .AND. &
-                        gp_outside_of_building(3) == 1 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz
-                      ei(num_gp)     = e(k+1,j,i)
-                      dissi(num_gp)  = diss(k+1,j,i)
-                      de_dxi(num_gp) = de_dx(k+1,j,i)
-                      de_dyi(num_gp) = de_dy(k+1,j,i)
-                      de_dzi(num_gp) = 0.0_wp
-                   ENDIF
-
-!
-!--                If wall between gridpoint 5 and gridpoint 7, then
-!--                Neumann boundary condition has to be applied
-!--                (only one case as only building beneath is possible)
-                   IF ( gp_outside_of_building(5) == 0  .AND. &
-                        gp_outside_of_building(7) == 1 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = j * dy
-                      location(num_gp,3) = k * dz
-                      ei(num_gp)     = e(k+1,j,i+1)
-                      dissi(num_gp)  = diss(k+1,j,i+1)
-                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
-                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
-                      de_dzi(num_gp) = 0.0_wp
-                   ENDIF
-
-!
-!--                If wall between gridpoint 2 and gridpoint 4, then
-!--                Neumann boundary condition has to be applied
-!--                (only one case as only building beneath is possible)
-                   IF ( gp_outside_of_building(2) == 0  .AND. &
-                        gp_outside_of_building(4) == 1 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = i * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz
-                      ei(num_gp)     = e(k+1,j+1,i)
-                      dissi(num_gp)  = diss(k+1,j+1,i)
-                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
-                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
-                      de_dzi(num_gp) = 0.0_wp
-                   ENDIF
-
-! 
-!--                If wall between gridpoint 6 and gridpoint 8, then
-!--                Neumann boundary condition has to be applied
-!--                (only one case as only building beneath is possible)
-                   IF ( gp_outside_of_building(6) == 0  .AND. &
-                        gp_outside_of_building(8) == 1 )  THEN
-                      num_gp = num_gp + 1
-                      location(num_gp,1) = (i+1) * dx
-                      location(num_gp,2) = (j+1) * dy
-                      location(num_gp,3) = k * dz
-                      ei(num_gp)     = e(k+1,j+1,i+1)
-                      dissi(num_gp)  = diss(k+1,j+1,i+1)
-                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
-                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
-                      de_dzi(num_gp) = 0.0_wp
-                   ENDIF
-     
-!
-!--                Carry out the interpolation
-                   IF ( num_gp == 1 )  THEN
-! 
-!--                   If only one of the gridpoints is situated outside of the
-!--                   building, it follows that the values at the particle
-!--                   location are the same as the gridpoint values
-                      e_int(n)     = ei(num_gp)    
-                      diss_int(n)  = dissi(num_gp) 
-                      de_dx_int(n) = de_dxi(num_gp)
-                      de_dy_int(n) = de_dyi(num_gp)
-                      de_dz_int(n) = de_dzi(num_gp)
-!
-!--                   Set flag for stochastic equation which disables calculation 
-!--                   of drift and memory term near topography.
-                      term_1_2(n) = 0.0_wp
-                   ELSE IF ( num_gp > 1 )  THEN
-    
-                      d_sum = 0.0_wp
-! 
-!--                   Evaluation of the distances between the gridpoints
-!--                   contributing to the interpolated values, and the particle
-!--                   location
-                      DO  agp = 1, num_gp
-                         d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
-                                      + ( particles(n)%y-location(agp,2) )**2  &
-                                      + ( zv(n)-location(agp,3) )**2
-                         d_sum        = d_sum + d_gp_pl(agp)
-                      ENDDO
-    
-!
-!--                  Finally the interpolation can be carried out
-                      e_int(n)     = 0.0_wp 
-                      diss_int(n)  = 0.0_wp 
-                      de_dx_int(n) = 0.0_wp  
-                      de_dy_int(n) = 0.0_wp  
-                      de_dz_int(n) = 0.0_wp  
-                      DO  agp = 1, num_gp
-                         e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
-                                                ei(agp) / ( (num_gp-1) * d_sum )
-                         diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
-                                             dissi(agp) / ( (num_gp-1) * d_sum )
-                         de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                            de_dxi(agp) / ( (num_gp-1) * d_sum )
-                         de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                            de_dyi(agp) / ( (num_gp-1) * d_sum )
-                         de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
-                                            de_dzi(agp) / ( (num_gp-1) * d_sum )
-                      ENDDO
-    
-                   ENDIF
-                   e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
-                   diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
-                   de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
-                   de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) ) 
-                   de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) ) 
-!
-!--                Set flag for stochastic equation which disables calculation 
-!--                of drift and memory term near topography.
-                   term_1_2(n) = 0.0_wp
-                ENDIF
-             ENDDO
-          ENDDO
-       ENDIF
+          ENDIF
+       ENDDO
 
        DO nb = 0,7
Index: palm/trunk/SOURCE/lpm_boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 2697)
+++ palm/trunk/SOURCE/lpm_boundary_conds.f90	(revision 2698)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Particle reflections at downward-facing walls implemented. Moreover, 
+! reflections are adjusted to revised particle grid box location. 
+! (responsible Philipp Thiele) 
 ! 
 ! Former revisions:
@@ -92,5 +94,5 @@
 !> (see offset_ocean_*)
 !------------------------------------------------------------------------------!
- SUBROUTINE lpm_boundary_conds( location )
+ SUBROUTINE lpm_boundary_conds( location , i, j, k )
  
 
@@ -108,5 +110,5 @@
 
     USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys, nz, nzb
+        ONLY:  nxl, nxr, nyn, nys, nz, nzb, wall_flags_0,nyng,nysg
 
     USE kinds
@@ -119,10 +121,11 @@
     USE pegrid
 
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index
-
     IMPLICIT NONE
 
     CHARACTER (LEN=*) ::  location     !<
+    
+    INTEGER(iwp), INTENT(IN) ::  i !< 
+    INTEGER(iwp), INTENT(IN) ::  j !<
+    INTEGER(iwp), INTENT(IN) ::  k !<
     
     INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
@@ -133,19 +136,23 @@
     INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
     INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
-    INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
-    INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
-    INTEGER(iwp) ::  k_wall         !< vertical index of topography top
+    INTEGER(iwp) ::  j2             !< grid index (y) of current particle position
+    INTEGER(iwp) ::  j3             !< grid index (y) of intermediate particle position
+    INTEGER(iwp) ::  k1             !< grid index (z) of old particle position
+    INTEGER(iwp) ::  k2             !< grid index (z) of current particle position
+    INTEGER(iwp) ::  k3             !< grid index (z) of intermediate particle position
     INTEGER(iwp) ::  n              !< particle number
     INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
     INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
-    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
-    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
+    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithm
+    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithm
+    INTEGER(iwp) ::  tmp_z          !< dummy for sorting algorithm
 
     INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions 
-    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions 
+    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
+    INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
     
     LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
     LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
-    LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
+    LOGICAL  ::  cross_wall_z    !< flag to check if particle reflection along z is necessary
     LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
     LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
@@ -156,4 +163,5 @@
     LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
     LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
+    LOGICAL  ::  z_wall_reached  !< flag to check if particle has already reached wall
 
     LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
@@ -177,8 +185,5 @@
     REAL(wp) ::  xwall          !< location of wall in x
     REAL(wp) ::  ywall          !< location of wall in y
-    REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
-    REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
-    REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
-    REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
+    REAL(wp) ::  zwall          !< location of wall in z
 
     REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
@@ -263,4 +268,6 @@
           i2 = particles(n)%x * ddx
           j2 = particles(n)%y * ddy
+          IF (zw(k)   < particles(n)%z ) k2 = k + 1
+          IF (zw(k-1) > particles(n)%z ) k2 = k - 1 
 !
 !--       Save current particle positions
@@ -275,6 +282,7 @@
 !
 !--       Obtain x/y indices for old particle positions
-          i1 = pos_x_old * ddx
-          j1 = pos_y_old * ddy
+          i1 = i
+          j1 = j
+          k1 = k
 !
 !--       Determine horizontal as well as vertical walls at which particle can 
@@ -283,9 +291,9 @@
 !--       Wall to the right 
           IF ( prt_x > pos_x_old )  THEN
-             xwall = ( i1 + 0.5_wp ) * dx
+             xwall = ( i1 + 1 ) * dx
 !
 !--       Wall to the left
           ELSE
-             xwall = ( i1 - 0.5_wp ) * dx
+             xwall = i1 * dx
           ENDIF
 !
@@ -293,22 +301,20 @@
 !--       Wall to the north
           IF ( prt_y > pos_y_old )  THEN
-             ywall = ( j1 + 0.5_wp ) * dy
+             ywall = ( j1 +1 ) * dy
 !--       Wall to the south
           ELSE
-             ywall = ( j1 - 0.5_wp ) * dy
-          ENDIF
-!
-!--       Walls aligned in xy layer at which particle can be possiblly reflected.
-!--       The construct of MERGE and BTEST is used to determine the topography-
-!--       top index (former nzb_s_inner). 
-          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' ) )
+             ywall = j1 * dy
+          ENDIF
+
+          IF ( prt_z > pos_z_old ) THEN
+             zwall = zw(k)
+          ELSE
+             zwall = zw(k-1)
+          ENDIF     
 !
 !--       Initialize flags to check if particle reflection is necessary
-          downwards    = .FALSE.
           cross_wall_x = .FALSE.
           cross_wall_y = .FALSE.
+          cross_wall_z = .FALSE.
 !
 !--       Initialize flags to check if a wall is reached
@@ -318,12 +324,13 @@
 !
 !--       Initialize flags to check if a particle was already reflected
-          reflect_x = .FALSE.
-          reflect_y = .FALSE.
-          reflect_z = .FALSE.
-!
-!--       Initialize flags to check if a vertical wall is already crossed.
+          reflect_x    = .FALSE.
+          reflect_y    = .FALSE.
+          reflect_z    = .FALSE.
+!
+!--       Initialize flags to check if a wall is already crossed.
 !--       ( Required to obtain correct indices. )
           x_wall_reached = .FALSE.
           y_wall_reached = .FALSE.
+          z_wall_reached = .FALSE.
 !
 !--       Initialize time array 
@@ -342,4 +349,5 @@
           x_ind(t_index)   = i2
           y_ind(t_index)   = j1
+          z_ind(t_index)   = k1
           reach_x(t_index) = .TRUE.
           reach_y(t_index) = .FALSE.
@@ -360,4 +368,5 @@
           x_ind(t_index)   = i1
           y_ind(t_index)   = j2
+          z_ind(t_index)   = k1
           reach_x(t_index) = .FALSE.
           reach_y(t_index) = .TRUE.
@@ -369,54 +378,24 @@
 !
 !--       z-direction
-!--       At first, check if particle moves downwards. Only in this case a 
-!--       particle can be reflected vertically.
-          IF ( prt_z < pos_z_old )  THEN
-
-             downwards = .TRUE.
-             dum       =  1.0_wp / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ), &
-                                          MIN( prt_z - pos_z_old, -1E-30_wp ), &
-                                          prt_z > pos_z_old )
-
-             t(t_index)       = ( zwall1 - pos_z_old ) * dum
-             x_ind(t_index)   = i2
-             y_ind(t_index)   = j2
-             reach_x(t_index) = .FALSE.
-             reach_y(t_index) = .FALSE.
-             reach_z(t_index) = .TRUE.
-             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
-                t_index = t_index + 1
-
-             reach_x(t_index) = .FALSE.
-             reach_y(t_index) = .FALSE.
-             reach_z(t_index) = .TRUE.
-             t(t_index)       = ( zwall2 - pos_z_old ) * dum
-             x_ind(t_index)   = i1
-             y_ind(t_index)   = j1
-             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
-                t_index = t_index + 1
-
-             reach_x(t_index) = .FALSE.
-             reach_y(t_index) = .FALSE.
-             reach_z(t_index) = .TRUE.
-             t(t_index)       = ( zwall3 - pos_z_old ) * dum
-             x_ind(t_index)   = i2
-             y_ind(t_index)   = j1
-             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
-                t_index = t_index + 1
-
-             reach_x(t_index) = .FALSE.
-             reach_y(t_index) = .FALSE.
-             reach_z(t_index) = .TRUE.
-             t(t_index)       = ( zwall4 - pos_z_old ) * dum
-             x_ind(t_index)   = i1
-             y_ind(t_index)   = j2
-             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
-                t_index = t_index + 1
-
-          END IF
+          t(t_index) = (zwall - pos_z_old )                                    &
+                     / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ),             &
+                              MIN( prt_z - pos_z_old, -1E-30_wp ),             &
+                              prt_z > pos_z_old )
+                     
+          x_ind(t_index)   = i1
+          y_ind(t_index)   = j1
+          z_ind(t_index)   = k2
+          reach_x(t_index) = .FALSE.
+          reach_y(t_index) = .FALSE.
+          reach_z(t_index) = .TRUE.
+          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN
+             t_index      = t_index + 1
+             cross_wall_z = .TRUE.
+          ENDIF
+          
           t_index_number = t_index - 1
 !
 !--       Carry out reflection only if particle reaches any wall
-          IF ( cross_wall_x .OR. cross_wall_y .OR. downwards )  THEN
+          IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z )  THEN
 !
 !--          Sort fractional timesteps in ascending order. Also sort the 
@@ -435,4 +414,5 @@
                    tmp_x       = x_ind(ir)
                    tmp_y       = y_ind(ir)
+                   tmp_z       = z_ind(ir)
                    tmp_reach_x = reach_x(ir)
                    tmp_reach_y = reach_y(ir)
@@ -443,4 +423,5 @@
                       x_ind(jr)   = x_ind(jr-inc)
                       y_ind(jr)   = y_ind(jr-inc)
+                      z_ind(jr)   = z_ind(jr-inc)
                       reach_x(jr) = reach_x(jr-inc)
                       reach_y(jr) = reach_y(jr-inc)
@@ -452,4 +433,5 @@
                    x_ind(jr)   = tmp_x
                    y_ind(jr)   = tmp_y
+                   z_ind(jr)   = tmp_z
                    reach_x(jr) = tmp_reach_x
                    reach_y(jr) = tmp_reach_y
@@ -480,24 +462,26 @@
                 i3 = x_ind(t_index)
                 j3 = y_ind(t_index)
+                k3 = z_ind(t_index)
 !
 !--             Check which wall is already reached
                 IF ( .NOT. x_wall_reached )  x_wall_reached = reach_x(t_index) 
-                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index) 
+                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index)
+                IF ( .NOT. z_wall_reached )  z_wall_reached = reach_z(t_index)
 !
 !--             Check if a particle needs to be reflected at any yz-wall. If 
 !--             necessary, carry out reflection. Please note, a security 
-!--             constant is required, as the particle position do not 
+!--             constant is required, as the particle position does not 
 !--             necessarily exactly match the wall location due to rounding 
-!--             errors. At first, determine index of topography top at (j3,i3)  
-                k_wall = get_topography_top_index( j3, i3, 's' ) 
-                IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
-                     pos_z <= zw(k_wall)             .AND.                     &
-                     reach_x(t_index)                .AND.                     &
+!--             errors.
+                IF ( reach_x(t_index)                      .AND.               & 
+                     ABS( pos_x - xwall ) < eps            .AND.               &
+                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
                      .NOT. reflect_x )  THEN
-!
+! 
+! 
 !--                Reflection in x-direction. 
 !--                Ensure correct reflection by MIN/MAX functions, depending on
 !--                direction of particle transport. 
-!--                Due to rounding errors pos_x do not exactly matches the wall
+!--                Due to rounding errors pos_x does not exactly match the wall
 !--                location, leading to erroneous reflection.             
                    pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
@@ -508,5 +492,5 @@
                    particles(n)%speed_x = - particles(n)%speed_x
 !
-!--                Change also sign of subgrid-scale particle speed
+!--                Also change sign of subgrid-scale particle speed
                    particles(n)%rvar1 = - particles(n)%rvar1
 !
@@ -514,5 +498,5 @@
                    reflect_x          = .TRUE.
 !
-!--                As particle do not crosses any further yz-wall during 
+!--                As the particle does not cross any further yz-wall during 
 !--                this timestep, set further x-indices to the current one.
                    x_ind(t_index:t_index_number) = i1
@@ -522,50 +506,90 @@
                 ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
                     x_ind(t_index:t_index_number) = i2
-                ENDIF
+                ENDIF !particle reflection in x direction done
+
 !
 !--             Check if a particle needs to be reflected at any xz-wall. If 
-!--             necessary, carry out reflection. At first, determine index of
-!--             topography top at (j3,i3)  
-                k_wall = get_topography_top_index( j3, i3, 's' ) 
-                IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
-                     pos_z <= zw(k_wall)             .AND.                     &
-                     reach_y(t_index)                .AND.                     &
-                     .NOT. reflect_y ) THEN
-
+!--             necessary, carry out reflection. Please note, a security 
+!--             constant is required, as the particle position does not 
+!--             necessarily exactly match the wall location due to rounding 
+!--             errors. 
+                WRITE(9,*) 'wall_test_y nysg:',nysg ,' y: ',particles(n)%y,' j3: ',j3,'nyng:', nyng  
+                FLUSH(9)
+                WRITE(9,*) BTEST(wall_flags_0(k3,j3,i3),0) 
+                FLUSH(9)
+                IF ( reach_y(t_index)                      .AND.               & 
+                     ABS( pos_y - ywall ) < eps            .AND.               &
+                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
+                     .NOT. reflect_y )  THEN
+! 
+! 
+!--                Reflection in y-direction. 
+!--                Ensure correct reflection by MIN/MAX functions, depending on
+!--                direction of particle transport. 
+!--                Due to rounding errors pos_y does not exactly match the wall
+!--                location, leading to erroneous reflection.             
                    pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
                                   MAX( 2.0_wp * ywall - pos_y, ywall ),        &
-                                  particles(n)%y > ywall ) 
-
-                   particles(n)%speed_y = - particles(n)%speed_y      
-                   particles(n)%rvar2   = - particles(n)%rvar2       
-     
-                   reflect_y            = .TRUE.
+                                  particles(n)%y > ywall )
+!
+!--                Change sign of particle speed                     
+                   particles(n)%speed_y = - particles(n)%speed_y
+!
+!--                Also change sign of subgrid-scale particle speed
+                   particles(n)%rvar2 = - particles(n)%rvar2
+!
+!--                Set flag that reflection along y is already done
+                   reflect_y          = .TRUE.
+!
+!--                As the particle does not cross any further xz-wall during 
+!--                this timestep, set further y-indices to the current one.
                    y_ind(t_index:t_index_number) = j1
-
+!
+!--             If particle already reached the wall but was not reflected, 
+!--             set further y-indices to the new one.
                 ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
-                   y_ind(t_index:t_index_number) = j2
-                ENDIF
+                    y_ind(t_index:t_index_number) = j2
+                ENDIF !particle reflection in y direction done
+                
 !
 !--             Check if a particle needs to be reflected at any xy-wall. If 
-!--             necessary, carry out reflection. 
-                IF ( downwards .AND. reach_z(t_index) .AND.                    &
+!--             necessary, carry out reflection. Please note, a security 
+!--             constant is required, as the particle position does not 
+!--             necessarily exactly match the wall location due to rounding 
+!--             errors. 
+                IF ( reach_z(t_index)                      .AND.               & 
+                     ABS( pos_z - zwall ) < eps            .AND.               &
+                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
                      .NOT. reflect_z )  THEN
-!
-!--                Determine index of topography top at (j3,i3) and chick if
-!--                particle is below.  
-                   k_wall = get_topography_top_index( j3, i3, 's' ) 
-                   IF ( pos_z - zw(k_wall) < eps ) THEN
-  
-                      pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z,    &
-                                   zw(k_wall) )
-
-                      particles(n)%speed_z = - particles(n)%speed_z
-                      particles(n)%rvar3   = - particles(n)%rvar3
-
-                      reflect_z            = .TRUE.
-
-                   ENDIF
-
-                ENDIF
+! 
+! 
+!--                Reflection in z-direction. 
+!--                Ensure correct reflection by MIN/MAX functions, depending on
+!--                direction of particle transport. 
+!--                Due to rounding errors pos_z does not exactly match the wall
+!--                location, leading to erroneous reflection.             
+                   pos_z = MERGE( MIN( 2.0_wp * zwall - pos_z, zwall ),        &
+                                  MAX( 2.0_wp * zwall - pos_z, zwall ),        &
+                                  particles(n)%z > zwall )
+!
+!--                Change sign of particle speed                     
+                   particles(n)%speed_z = - particles(n)%speed_z
+!
+!--                Also change sign of subgrid-scale particle speed
+                   particles(n)%rvar3 = - particles(n)%rvar3
+!
+!--                Set flag that reflection along z is already done
+                   reflect_z          = .TRUE.
+!
+!--                As the particle does not cross any further xy-wall during 
+!--                this timestep, set further z-indices to the current one.
+                   z_ind(t_index:t_index_number) = k1
+!
+!--             If particle already reached the wall but was not reflected, 
+!--             set further z-indices to the new one.
+                ELSEIF ( z_wall_reached .AND. .NOT. reflect_z )  THEN
+                    z_ind(t_index:t_index_number) = k2
+                ENDIF !particle reflection in z direction done                
+                
 !
 !--             Swap time
Index: palm/trunk/SOURCE/lpm_init.f90
===================================================================
--- palm/trunk/SOURCE/lpm_init.f90	(revision 2697)
+++ palm/trunk/SOURCE/lpm_init.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele)
 ! 
 ! Former revisions:
@@ -236,5 +236,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -757,5 +757,5 @@
 !--                            Determine surface level. Therefore, check for
 !--                            upward-facing wall on w-grid. 
-                               k_surf = get_topography_top_index( jp, ip, 'w' )
+                               k_surf = get_topography_top_index_ji( jp, ip, 'w' )
 
                                IF ( seed_follows_topography )  THEN
@@ -930,5 +930,5 @@
 !--             Identify particles located outside the model domain and reflect
 !--             or absorb them if necessary.
-                CALL lpm_boundary_conds( 'bottom/top' )
+                CALL lpm_boundary_conds( 'bottom/top', i, j, k )
 !
 !--             Furthermore, remove particles located in topography. Note, as
Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -928,5 +928,5 @@
 
        USE surface_mod,                                                        &
-           ONLY:  get_topography_top_index
+           ONLY:  get_topography_top_index_ji
 
 
@@ -945,5 +945,5 @@
 !
 !--          Determine vertical index of topography top
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
              DO  k = nzb+1, nzt
 !
@@ -2319,5 +2319,5 @@
 
        USE surface_mod,                                                        &
-           ONLY:  get_topography_top_index
+           ONLY:  get_topography_top_index_ji
 
 
@@ -2334,5 +2334,5 @@
 !
 !--    Determine vertical index of topography top
-       k_wall = get_topography_top_index( j, i, 's' )
+       k_wall = get_topography_top_index_ji( 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 2697)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -151,5 +151,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index
+        ONLY:  get_topography_top_index_ji
 
 
@@ -400,5 +400,5 @@
                DO  j = nys, nyn
                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index( j, i, 's' )
+                     k_topo = get_topography_top_index_ji( j, i, 's' )
                      DO  k = k_topo, k_topo + pch_index_ji(j,i)
                         local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
@@ -415,5 +415,5 @@
                DO  j = nys, nyn
                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index( j, i, 's' )
+                     k_topo = get_topography_top_index_ji( j, i, 's' )
                      DO  k = k_topo, k_topo + pch_index_ji(j,i)
                         local_pf(i,j,k) = lad_s(k-k_topo,j,i)
@@ -803,5 +803,5 @@
 !--          Check whether topography and local vegetation on top exceed 
 !--          height of the model domain.
-             k = get_topography_top_index( j, i, 's' )
+             k = get_topography_top_index_ji( j, i, 's' )
              IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
                 message_string =  'Local vegetation height on top of ' //      &
@@ -1110,5 +1110,5 @@
 !
 !--                Determine topography-top index on u-grid
-                   k_wall = get_topography_top_index( j, i, 'u' )
+                   k_wall = get_topography_top_index_ji( j, i, 'u' )
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
@@ -1174,5 +1174,5 @@
 !
 !--                Determine topography-top index on v-grid
-                   k_wall = get_topography_top_index( j, i, 'v' )
+                   k_wall = get_topography_top_index_ji( j, i, 'v' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1239,5 +1239,5 @@
 !
 !--                Determine topography-top index on w-grid
-                   k_wall = get_topography_top_index( j, i, 'w' )
+                   k_wall = get_topography_top_index_ji( j, i, 'w' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
@@ -1291,5 +1291,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index( j, i, 's' )
+                   k_wall = get_topography_top_index_ji( j, i, 's' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1308,5 +1308,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index( j, i, 's' )
+                   k_wall = get_topography_top_index_ji( j, i, 's' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1338,5 +1338,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index( j, i, 's' )
+                   k_wall = get_topography_top_index_ji( j, i, 's' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1367,5 +1367,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index( j, i, 's' )
+                   k_wall = get_topography_top_index_ji( j, i, 's' )
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1461,5 +1461,5 @@
 !
 !--          Determine topography-top index on u-grid
-             k_wall = get_topography_top_index( j, i, 'u' )
+             k_wall = get_topography_top_index_ji( j, i, 'u' )
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
@@ -1521,5 +1521,5 @@
 !
 !--          Determine topography-top index on v-grid
-             k_wall = get_topography_top_index( j, i, 'v' )
+             k_wall = get_topography_top_index_ji( j, i, 'v' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -1581,5 +1581,5 @@
 !
 !--          Determine topography-top index on w-grid
-             k_wall = get_topography_top_index( j, i, 'w' )
+             k_wall = get_topography_top_index_ji( j, i, 'w' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
@@ -1628,5 +1628,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -1641,5 +1641,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -1667,5 +1667,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -1693,5 +1693,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -265,5 +265,5 @@
 
     USE surface_mod,                                                            &
-        ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -1410,17 +1410,17 @@
 !--             Determine largest topography index on scalar grid
                 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
-                                         get_topography_top_index( j, i, 's' ) )
+                                         get_topography_top_index_ji( j, i, 's' ) )
 !
 !--             Determine largest topography index on u grid
                 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
-                                         get_topography_top_index( j, i, 'u' ) )
+                                         get_topography_top_index_ji( j, i, 'u' ) )
 !
 !--             Determine largest topography index on v grid
                 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
-                                         get_topography_top_index( j, i, 'v' ) )
+                                         get_topography_top_index_ji( j, i, 'v' ) )
 !
 !--             Determine largest topography index on w grid
                 nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
-                                         get_topography_top_index( j, i, 'w' ) )
+                                         get_topography_top_index_ji( j, i, 'w' ) )
              ENDDO
           ENDDO
@@ -1436,17 +1436,17 @@
 !--             Determine largest topography index on scalar grid
                 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
-                                         get_topography_top_index( j, i, 's' ) )
+                                         get_topography_top_index_ji( j, i, 's' ) )
 !
 !--             Determine largest topography index on u grid
                 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
-                                         get_topography_top_index( j, i, 'u' ) )
+                                         get_topography_top_index_ji( j, i, 'u' ) )
 !
 !--             Determine largest topography index on v grid
                 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
-                                         get_topography_top_index( j, i, 'v' ) )
+                                         get_topography_top_index_ji( j, i, 'v' ) )
 !
 !--             Determine largest topography index on w grid
                 nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
-                                         get_topography_top_index( j, i, 'w' ) )
+                                         get_topography_top_index_ji( j, i, 'w' ) )
           ENDDO
           nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
@@ -1461,17 +1461,17 @@
 !--             Determine largest topography index on scalar grid
                 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
-                                         get_topography_top_index( j, i, 's' ) )
+                                         get_topography_top_index_ji( j, i, 's' ) )
 !
 !--             Determine largest topography index on u grid
                 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
-                                         get_topography_top_index( j, i, 'u' ) )
+                                         get_topography_top_index_ji( j, i, 'u' ) )
 !
 !--             Determine largest topography index on v grid
                 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
-                                         get_topography_top_index( j, i, 'v' ) )
+                                         get_topography_top_index_ji( j, i, 'v' ) )
 !
 !--             Determine largest topography index on w grid
                 nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
-                                         get_topography_top_index( j, i, 'w' ) )
+                                         get_topography_top_index_ji( j, i, 'w' ) )
              ENDDO
           ENDDO
@@ -1487,17 +1487,17 @@
 !--             Determine largest topography index on scalar grid
                 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
-                                         get_topography_top_index( j, i, 's' ) )
+                                         get_topography_top_index_ji( j, i, 's' ) )
 !
 !--             Determine largest topography index on u grid
                 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
-                                         get_topography_top_index( j, i, 'u' ) )
+                                         get_topography_top_index_ji( j, i, 'u' ) )
 !
 !--             Determine largest topography index on v grid
                 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
-                                         get_topography_top_index( j, i, 'v' ) )
+                                         get_topography_top_index_ji( j, i, 'v' ) )
 !
 !--             Determine largest topography index on w grid
                 nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
-                                         get_topography_top_index( j, i, 'w' ) )
+                                         get_topography_top_index_ji( j, i, 'w' ) )
           ENDDO
           nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
@@ -1552,5 +1552,5 @@
 !--          is part of the surfacetypes now. Set default roughness instead. 
 !--          Determine topography top index on u-grid
-             kb  = get_topography_top_index( j, i, 'u' )
+             kb  = get_topography_top_index_ji( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1568,5 +1568,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = get_topography_top_index( j, i, 'v' )
+             kb  = get_topography_top_index_ji( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1611,5 +1611,5 @@
 !--          to the present surface tpye. 
 !--          Determine topography top index on u-grid
-             kb  = get_topography_top_index( j, i, 'u' )
+             kb  = get_topography_top_index_ji( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1627,5 +1627,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = get_topography_top_index( j, i, 'v' )
+             kb  = get_topography_top_index_ji( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1667,5 +1667,5 @@
 !
 !--          Determine topography top index on u-grid
-             kb  = get_topography_top_index( j, i, 'u' )
+             kb  = get_topography_top_index_ji( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1683,5 +1683,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = get_topography_top_index( j, i, 'v' )
+             kb  = get_topography_top_index_ji( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1723,5 +1723,5 @@
 !
 !--          Determine topography top index on u-grid
-             kb  = get_topography_top_index( j, i, 'u' )
+             kb  = get_topography_top_index_ji( j, i, 'u' )
              k   = kb + 1
              wall_index = kb
@@ -1739,5 +1739,5 @@
 !
 !--          Determine topography top index on v-grid
-             kb  = get_topography_top_index( j, i, 'v' )
+             kb  = get_topography_top_index_ji( j, i, 'v' )
              k   = kb + 1
              wall_index = kb
@@ -1779,11 +1779,11 @@
 !
 !--             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' )
+                k_wall_u_ji   = get_topography_top_index_ji( j,   0, 'u_out' )
+                k_wall_u_ji_p = get_topography_top_index_ji( j+1, 0, 'u_out' )
+                k_wall_u_ji_m = get_topography_top_index_ji( j-1, 0, 'u_out' )
+
+                k_wall_w_ji   = get_topography_top_index_ji( j,   -1, 'w_out' )
+                k_wall_w_ji_p = get_topography_top_index_ji( j+1, -1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index_ji( j-1, -1, 'w_out' )
 
                 DO  k = nzb, nzt_topo_nestbc_l
@@ -1869,11 +1869,11 @@
 !
 !--             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' )
+                k_wall_u_ji   = get_topography_top_index_ji( j,   i, 'u_out' )
+                k_wall_u_ji_p = get_topography_top_index_ji( j+1, i, 'u_out' )
+                k_wall_u_ji_m = get_topography_top_index_ji( j-1, i, 'u_out' )
+
+                k_wall_w_ji   = get_topography_top_index_ji( j,   i, 'w_out' )
+                k_wall_w_ji_p = get_topography_top_index_ji( j+1, i, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index_ji( j-1, i, 'w_out' )
 
                 DO  k = nzb, nzt_topo_nestbc_r
@@ -1952,11 +1952,11 @@
 !
 !--             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' )
+                k_wall_v_ji   = get_topography_top_index_ji( 0, i,   'v_out' )
+                k_wall_v_ji_p = get_topography_top_index_ji( 0, i+1, 'v_out' )
+                k_wall_v_ji_m = get_topography_top_index_ji( 0, i-1, 'v_out' )
+
+                k_wall_w_ji   = get_topography_top_index_ji( -1, i,   'w_out' )
+                k_wall_w_ji_p = get_topography_top_index_ji( -1, i+1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index_ji( -1, i-1, 'w_out' )
 
                 DO  k = nzb, nzt_topo_nestbc_s
@@ -2041,11 +2041,11 @@
 !
 !--             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' )
+                k_wall_v_ji   = get_topography_top_index_ji( j, i,   'v_out' )
+                k_wall_v_ji_p = get_topography_top_index_ji( j, i+1, 'v_out' )
+                k_wall_v_ji_m = get_topography_top_index_ji( j, i-1, 'v_out' )
+
+                k_wall_w_ji   = get_topography_top_index_ji( j, i,   'w_out' )
+                k_wall_w_ji_p = get_topography_top_index_ji( j, i+1, 'w_out' )
+                k_wall_w_ji_m = get_topography_top_index_ji( j, i-1, 'w_out' )
 
                 DO  k = nzb, nzt_topo_nestbc_n
@@ -2643,5 +2643,5 @@
           i = nxl - 1
           DO  j = nysg, nyng
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2664,5 +2664,5 @@
           i = nxr + 1
           DO  j = nysg, nyng
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2685,5 +2685,5 @@
           j = nys - 1
           DO  i = nxlg, nxrg
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
              
              DO  k = k_wall + 1, nzt
@@ -2706,5 +2706,5 @@
           j = nyn + 1
           DO  i = nxlg, nxrg
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              DO  k = k_wall + 1, nzt
@@ -2729,5 +2729,5 @@
 !
 !--          Determine vertical index for local topography top
-             k_wall = get_topography_top_index( j, i, 's' )
+             k_wall = get_topography_top_index_ji( j, i, 's' )
 
              kc     = kco(k) + 1
@@ -3196,5 +3196,5 @@
 !
 !--             Determine vertical index of topography top at grid point (j,i)
-                k_wall = get_topography_top_index( j, i, TRIM ( var ) )
+                k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
 !
 !--             kbc is the first coarse-grid point above the surface
@@ -3227,5 +3227,5 @@
 !
 !--              Determine vertical index of topography top at grid point (j,i)
-                 k_wall = get_topography_top_index( j, i, 'w' )
+                 k_wall = get_topography_top_index_ji( j, i, 'w' )
 
                  f(k_wall,j,i) = 0.0_wp
@@ -4369,5 +4369,5 @@
 !
 !--      Determine vertical index of topography top at grid point (j,i)
-         k_wall = get_topography_top_index( j, i, TRIM( var ) )
+         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
          DO  k = k_wall, nzt+1
@@ -4395,5 +4395,5 @@
 !
 !--         Determine vertical index of topography top at grid point (j,i)
-            k_wall = get_topography_top_index( j, i, TRIM ( var ) )
+            k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
 
             k = k_wall+1
@@ -4420,5 +4420,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, TRIM ( var ) )
+               k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
 
                DO  k = k_wall+1, nzt_topo_nestbc
@@ -4445,5 +4445,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, TRIM ( var ) )
+               k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
 
                k = k_wall + 1
@@ -4476,5 +4476,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, 's' )
+               k_wall = get_topography_top_index_ji( j, i, 's' )
 
                DO  k = k_wall, nzt + 1
@@ -4486,5 +4486,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, 's' )
+               k_wall = get_topography_top_index_ji( j, i, 's' )
 
                DO  k = k_wall, nzt+1
@@ -4593,5 +4593,5 @@
 !
 !--      Determine vertical index of topography top at grid point (j,i)
-         k_wall = get_topography_top_index( j, i, TRIM( var ) )
+         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
          DO  k = k_wall, nzt+1
@@ -4619,5 +4619,5 @@
 !
 !--         Determine vertical index of topography top at grid point (j,i)
-            k_wall = get_topography_top_index( j, i, TRIM( var ) )
+            k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
             k = k_wall + 1
@@ -4643,5 +4643,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, TRIM( var ) )
+               k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
                DO  k = k_wall, nzt_topo_nestbc
@@ -4670,5 +4670,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, TRIM( var ) )
+               k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
                k = k_wall + 1
@@ -4701,5 +4701,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, 's' )
+               k_wall = get_topography_top_index_ji( j, i, 's' )
                DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
@@ -4710,5 +4710,5 @@
 !
 !--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index( j, i, 's' )
+               k_wall = get_topography_top_index_ji( j, i, 's' )
                DO  k = k_wall, nzt+1
                   f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
@@ -4891,5 +4891,5 @@
 !
 !--       Determine vertical index of topography top at grid point (j,i)
-          k_wall = get_topography_top_index( j, i, TRIM( var ) )
+          k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
           DO  k = k_wall, nzt+1
@@ -4970,5 +4970,5 @@
 !
 !--       Determine vertical index of topography top at grid point (j,i)
-          k_wall = get_topography_top_index( j, i, TRIM( var ) )
+          k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
 
           DO  k = k_wall, nzt+1
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -272,5 +272,6 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, surf_def_h, surf_def_v, surf_lsm_h,   &
+        ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
+               surf_def_h, surf_def_v, surf_lsm_h,                             &
                surf_lsm_v, surf_type, surf_usm_h, surf_usm_v
 
@@ -3244,5 +3245,5 @@
 !
 !--             Obtain topography top index (lower bound of RRTMG)
-                k_topo = get_topography_top_index( j, i, 's' )
+                k_topo = get_topography_top_index_ji( j, i, 's' )
 
                 IF ( lw_radiation )  THEN
@@ -4411,8 +4412,4 @@
                   pc_heating_rate, lad_s, prototype_lad, usm_lad_rma       
        
-       USE surface_mod,                                                        &
-           ONLY:  get_topography_top_index, surf_lsm_h, surf_lsm_v, surf_usm_h,&
-                  surf_usm_v
-
        IMPLICIT NONE
 
@@ -4448,5 +4445,5 @@
 !
 !--                Find topography top index
-                   k_topo = get_topography_top_index( j, i, 's' )
+                   k_topo = get_topography_top_index_ji( j, i, 's' )
 
                    DO k = nzt+1, 0, -1
@@ -4553,6 +4550,6 @@
                 DO  j = ijdb(3,ids), ijdb(4,ids)
 
-                   k_topo  = get_topography_top_index( j, i, 's' )
-                   k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
+                   k_topo  = get_topography_top_index_ji( j, i, 's' )
+                   k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
 
 
@@ -4572,5 +4569,5 @@
              DO j = nys, nyn
 !--              Find topography top index
-                 k_topo = get_topography_top_index( j, i, 's' )
+                 k_topo = get_topography_top_index_ji( j, i, 's' )
                  k = nzut - k_topo
                  nsurfl = nsurfl + 6 * k
@@ -4593,5 +4590,5 @@
 !
 !--                Find topography top index
-                   k_topo = get_topography_top_index( j, i, 's' )
+                   k_topo = get_topography_top_index_ji( j, i, 's' )
 
                    DO k = k_topo + 1, pct(j,i)
@@ -4709,6 +4706,6 @@
                DO i = ijdb(1,ids), ijdb(2,ids)
                    DO j = ijdb(3,ids), ijdb(4,ids)
-                       k_topo  = get_topography_top_index( j, i, 's' )
-                       k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
+                       k_topo  = get_topography_top_index_ji( j, i, 's' )
+                       k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
 
                        DO k = MAX(k_topo,k_topo2)+1, nzut
@@ -4727,5 +4724,5 @@
           DO i = nxl, nxr
              DO j = nys, nyn
-                k_topo = get_topography_top_index( j, i, 's' )
+                k_topo = get_topography_top_index_ji( j, i, 's' )
 
 !--             add upward surface
@@ -4746,5 +4743,5 @@
           DO i = nxl, nxr
              DO j = nys, nyn
-                k_topo = get_topography_top_index( j, i, 's' )
+                k_topo = get_topography_top_index_ji( j, i, 's' )
 !--             north
                 IF ( j /= ny ) THEN
@@ -4752,5 +4749,5 @@
                    jr = min(max(j-jdir(ids),0),ny)
                    ir = min(max(i-idir(ids),0),nx)
-                   k_topo2 = get_topography_top_index( jr, ir, 's' )
+                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
                    DO k = MAX(k_topo,k_topo2)+1, nzut
                       isurf = isurf + 1
@@ -4763,5 +4760,5 @@
                    jr = min(max(j-jdir(ids),0),ny)
                    ir = min(max(i-idir(ids),0),nx)
-                   k_topo2 = get_topography_top_index( jr, ir, 's' )
+                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
 
                    DO k = MAX(k_topo,k_topo2)+1, nzut
@@ -4775,5 +4772,5 @@
                    jr = min(max(j-jdir(ids),0),ny)
                    ir = min(max(i-idir(ids),0),nx)
-                   k_topo2 = get_topography_top_index( jr, ir, 's' )
+                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
 
                    DO k = MAX(k_topo,k_topo2)+1, nzut
@@ -4787,5 +4784,5 @@
                    jr = min(max(j-jdir(ids),0),ny)
                    ir = min(max(i-idir(ids),0),nx)
-                   k_topo2 = get_topography_top_index( jr, ir, 's' )
+                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
 
                    DO k = MAX(k_topo,k_topo2)+1, nzut
@@ -5255,5 +5252,5 @@
 !
 !--             Following expression equals former kk = k - nzb_s_inner(j,i) 
-                kk = k - get_topography_top_index( j, i, 's' )  !- lad arrays are defined flat
+                kk = k - get_topography_top_index_ji( 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
@@ -5856,5 +5853,5 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
-                    k = get_topography_top_index( j, i, 's' )
+                    k = get_topography_top_index_ji( j, i, 's' )
 
                     usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 2698)
@@ -391,9 +391,4 @@
     PRIVATE
 
-    INTERFACE get_topography_top_index
-       MODULE PROCEDURE get_topography_top_index
-       MODULE PROCEDURE get_topography_top_index_ji
-    END  INTERFACE get_topography_top_index
-
     INTERFACE init_bc
        MODULE PROCEDURE init_bc
@@ -431,5 +426,6 @@
 !
 !-- Public subroutines and functions
-    PUBLIC get_topography_top_index, init_bc, init_surfaces,                   &
+    PUBLIC get_topography_top_index, get_topography_top_index_ji, init_bc,     &
+           init_surfaces,                                                      &
            init_surface_arrays, surface_read_restart_data,                     & 
            surface_restore_elements, surface_write_restart_data,               &
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2697)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 2698)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in get_topography_top_index
 ! 
 ! Former revisions:
@@ -870,5 +870,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index
+        ONLY:  get_topography_top_index_ji
 
     IMPLICIT NONE
@@ -972,5 +972,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                nz_s_shift = get_topography_top_index( j, i, 's' )
+                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
 
                 e(nz_s_shift:nzt+1,j,i)  =  e(0:nzt+1-nz_s_shift,j,i)
@@ -991,5 +991,5 @@
           IF ( complex_terrain )  THEN
              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
-                nz_s_shift_l = get_topography_top_index( 0, 0, 's' )
+                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
              ELSE
                 nz_s_shift_l = 0
