Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4167)
+++ palm/trunk/SOURCE/Makefile	(revision 4168)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# Remove some dependencies on surface_mod that are no longer required without
+# get_topography_top_index functions
+# 
+# 4167 2019-08-16 11:01:48Z suehring
 # Remove no longer needed dependencies on surface_mod
 #
@@ -1237,5 +1241,4 @@
 	modules.o \
 	netcdf_interface_mod.o \
-	surface_mod.o \
 	random_function_mod.o
 mod_kinds.o: \
@@ -1272,6 +1275,5 @@
 	mod_kinds.o \
 	modules.o \
-	netcdf_data_input_mod.o \
-	surface_mod.o
+	netcdf_data_input_mod.o
 ocean_mod.o: \
 	advec_s_pw.o \
Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/biometeorology_mod.f90	(revision 4168)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4144 2019-08-06 09:11:47Z raasch
 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
 ! 
@@ -178,5 +181,5 @@
     USE indices,                                                               &
         ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
-               nysg, nyng
+               nysg, nyng, topo_top_ind
 
     USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
@@ -192,7 +195,4 @@
                radiation_interactions, rad_sw_in,                              &
                rad_sw_out, rad_lw_in, rad_lw_out
-
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
 
     IMPLICIT NONE
@@ -1557,6 +1557,5 @@
              j = mrtbl(iy,l)
              k = mrtbl(iz,l)
-             IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
-                   bio_cell_level + 1_iwp)  THEN
+             IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
 !
 !-- Averaging was done before, so we can just copy the result here
@@ -1583,6 +1582,5 @@
           j = mrtbl(iy,l)
           k = mrtbl(iz,l)
-          IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
-                bio_cell_level + 1_iwp)  THEN
+          IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
              IF ( mrt_include_sw )  THEN
                  tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
@@ -1608,5 +1606,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
-    pair, tmrt )
+                                            pair, tmrt )
 
     IMPLICIT NONE
@@ -1633,5 +1631,5 @@
 !-- Determine cell level closest to 1.1m above ground
 !   by making use of truncation due to int cast
-    k = INT( get_topography_top_index_ji(j, i, 's') + bio_cell_level )  !< Vertical cell center closest to 1.1m
+    k = INT( topo_top_ind(j,i,0) + bio_cell_level )  !< Vertical cell center closest to 1.1m
 
 !
Index: palm/trunk/SOURCE/bulk_cloud_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/bulk_cloud_model_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4110 2019-07-22 17:05:21Z suehring
 ! Pass integer flag array as well as boundary flags to WS scalar advection 
 ! routine
@@ -280,4 +283,5 @@
         ONLY:  advc_flags_s,                                                   &
                nbgp, nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt,     &
+               topo_top_ind,                                                   &
                wall_flags_0
 
@@ -291,5 +295,6 @@
 
     USE surface_mod,                                                           &
-        ONLY :  bc_h, get_topography_top_index_ji, surf_bulk_cloud_model,      &
+        ONLY :  bc_h,                                                          &
+                surf_bulk_cloud_model,                                         &
                 surf_microphysics_morrison, surf_microphysics_seifert, &
                 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
@@ -3807,5 +3812,5 @@
 !
 !--          Determine vertical index of topography top
-             k_wall = get_topography_top_index_ji( j, i, 's' )
+             k_wall = topo_top_ind(j,i,0)
              DO  k = nzb+1, nzt
 !
@@ -3854,5 +3859,5 @@
 !
 !--    Determine vertical index of topography top
-       k_wall = get_topography_top_index_ji( j, i, 's' )
+       k_wall = topo_top_ind(j,i,0)
        DO  k = nzb+1, nzt
 !
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 4167)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Remove variable grid
+! 
+! 4167 2019-08-16 11:01:48Z suehring
 ! Changed behaviour of masked output over surface to follow terrain and ignore 
 ! buildings (J.Resler, T.Gronemeier)
@@ -225,6 +228,4 @@
 
     IMPLICIT NONE
-
-    CHARACTER(LEN=5) ::  grid !< flag to distinquish between staggered grids
 
     INTEGER(iwp) ::  av                      !< flag for (non-)average output
@@ -310,7 +311,4 @@
        ENDIF
 !
-!--    Set default grid for terrain-following output
-       grid = 's'
-!
 !--    Store the variable chosen.
        resorted = .FALSE.
@@ -373,5 +371,6 @@
                          im = mask_i(mid,i)
                          jm = mask_j(mid,j)
-                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )),&
+                                       DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
                             kk = MIN( ktt+mask_k(mid,k), nzt+1 )
@@ -685,5 +684,4 @@
 
           CASE ( 'w' )
-             grid = 'w'
              IF ( av == 0 )  THEN
                 to_be_resorted => w
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 4167)
+++ palm/trunk/SOURCE/header.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration 
@@ -449,5 +452,5 @@
     USE indices,                                                               &
         ONLY:  mg_loc_ind, nnx, nny, nnz, nx, ny, nxl_mg, nxr_mg, nyn_mg,      &
-               nys_mg, nzt, nzt_mg
+               nys_mg, nzt, nzt_mg, topo_top_ind
 
     USE kinds
@@ -478,5 +481,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  surf_def_h, get_topography_top_index_ji
+        ONLY:  surf_def_h
 
     USE turbulence_closure_mod,                                                &
@@ -1111,5 +1114,5 @@
        WRITE( io, 280 ) 
        IF ( turbulent_inflow )  THEN
-          WRITE( io, 281 )  zu( get_topography_top_index_ji( 0, 0, 's' ) )
+          WRITE( io, 281 )  zu(topo_top_ind(0,0,0))
        ENDIF
        IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 4167)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4151 2019-08-09 08:24:30Z suehring
 ! Add netcdf directive around input calls (fix for last commit)
 ! 
@@ -719,6 +722,5 @@
                 surf_def_v,                                                    &
                 surf_lsm_h,                                                    &
-                surf_usm_h,                                                    &
-                get_topography_top_index_ji
+                surf_usm_h
 
 #if defined( _OPENACC )
@@ -1633,8 +1635,8 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                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' )
+                nz_u_shift = topo_top_ind(j,i,1)
+                nz_v_shift = topo_top_ind(j,i,2)
+                nz_w_shift = topo_top_ind(j,i,3)
+                nz_s_shift = topo_top_ind(j,i,0)
 
                 u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)                
@@ -1677,8 +1679,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_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' )
+                nz_u_shift_l = topo_top_ind(j,i,1)
+                nz_v_shift_l = topo_top_ind(j,i,2)
+                nz_w_shift_l = topo_top_ind(j,i,3)
+                nz_s_shift_l = topo_top_ind(j,i,0)
              ELSE
                 nz_u_shift_l = 0
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 4167)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 4168)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Pre-calculate topography top index and store it on an array (replaces former
+! functions get_topography_top_index)
+! 
+! 4159 2019-08-15 13:31:35Z suehring
 ! Revision of topography processing. This was not consistent between 2D and 3D
 ! buildings. 
@@ -420,5 +424,5 @@
                nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,                 &
                nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner,             &
-               nzb_w_outer, nzt, topo_min_level
+               nzb_w_outer, nzt, topo_top_ind, topo_min_level
     
     USE kinds
@@ -429,5 +433,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, get_topography_top_index_ji, init_bc
+        ONLY:  init_bc
 
     USE vertical_nesting_mod,                                                  &
@@ -919,8 +923,8 @@
     topo_min_level   = 0
 #if defined( __parallel )
-    CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),             &
+    CALL MPI_ALLREDUCE( MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),             &
                         topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
 #else
-    topo_min_level = MINVAL( get_topography_top_index( 's' ) )
+    topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 #endif
 !
@@ -960,9 +964,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_ji( j, i, 's' ) )
+             zu_s_inner(i,j) = zu(topo_top_ind(j,i,0))
 !
 !--          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_ji( j, i, 's' ) )
+             zw_w_inner(i,j) = zw(topo_top_ind(j,i,3))
           ENDDO
        ENDDO
@@ -988,5 +992,5 @@
 !
 !-- Initialize 2D-index arrays. Note, these will be removed soon! 
-    nzb_local(nys:nyn,nxl:nxr) = get_topography_top_index( 's' )
+    nzb_local(nys:nyn,nxl:nxr) = topo_top_ind(nys:nyn,nxl:nxr,0)
     CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
 !
@@ -996,8 +1000,8 @@
     IF ( TRIM( topography ) /= 'flat' )  THEN
 #if defined( __parallel )
-       CALL MPI_ALLREDUCE( MAXVAL( get_topography_top_index( 's' ) ),          &
+       CALL MPI_ALLREDUCE( MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ),          &
                            nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )               
 #else
-       nzb_local_max = MAXVAL( get_topography_top_index( 's' ) )
+       nzb_local_max = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 #endif
        nzb_local_min = topo_min_level
@@ -2017,7 +2021,4 @@
     USE pegrid
 
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, get_topography_top_index_ji
-
     IMPLICIT NONE
 
@@ -2572,5 +2573,5 @@
     USE indices,                                                               &
         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
-               nzt, wall_flags_0
+               nzt, topo_top_ind, wall_flags_0
 
     USE kinds
@@ -2579,4 +2580,5 @@
 
     INTEGER(iwp) ::  i             !< index variable along x
+    INTEGER(iwp) ::  ibit          !< integer bit position of topgraphy masking array
     INTEGER(iwp) ::  j             !< index variable along y
     INTEGER(iwp) ::  k             !< index variable along z
@@ -2885,5 +2887,48 @@
        ENDIF      
     ENDIF
-
+!
+!-- Pre-calculate topography top indices (former get_topography_top_index 
+!-- function)
+    ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:4) )
+!
+!-- Uppermost topography index on scalar grid
+    ibit = 12
+    topo_top_ind(:,:,0) = MAXLOC(                                              &
+                                  MERGE( 1, 0,                                 &
+                                          BTEST( wall_flags_0(:,:,:), ibit )   &
+                                       ), DIM = 1                              &
+                                ) - 1 
+!
+!-- Uppermost topography index on u grid 
+    ibit = 14
+    topo_top_ind(:,:,1) = MAXLOC(                                              &
+                                  MERGE( 1, 0,                                 &
+                                          BTEST( wall_flags_0(:,:,:), ibit )   &
+                                       ), DIM = 1                              &
+                                ) - 1 
+!
+!-- Uppermost topography index on v grid 
+    ibit = 16
+    topo_top_ind(:,:,2) = MAXLOC(                                              &
+                                  MERGE( 1, 0,                                 &
+                                          BTEST( wall_flags_0(:,:,:), ibit )   &
+                                       ), DIM = 1                              &
+                                ) - 1 
+!
+!-- Uppermost topography index on w grid
+    ibit = 18
+    topo_top_ind(:,:,3) = MAXLOC(                                              &
+                                  MERGE( 1, 0,                                 &
+                                          BTEST( wall_flags_0(:,:,:), ibit )   &
+                                       ), DIM = 1                              &
+                                ) - 1 
+!
+!-- Uppermost topography index on scalar outer grid
+    ibit = 24
+    topo_top_ind(:,:,4) = MAXLOC(                                              &
+                                  MERGE( 1, 0,                                 &
+                                          BTEST( wall_flags_0(:,:,:), ibit )   &
+                                       ), DIM = 1                              &
+                                ) - 1                            
 
  END SUBROUTINE set_topo_flags
Index: palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4145 2019-08-06 09:55:22Z schwenkel
 ! Some reformatting
 ! 
@@ -227,5 +230,7 @@
     USE indices,                                                               &
         ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb,    &
-               nzb_max, nzt, wall_flags_0,nbgp, ngp_2dh_outer
+               nzb_max, nzt,nbgp, ngp_2dh_outer,                               &
+               topo_top_ind,                                                   &
+               wall_flags_0
 
     USE kinds
@@ -259,6 +264,8 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h,&
-               bc_h
+        ONLY:  bc_h,                                                           &
+               surf_def_h,                                                     &
+               surf_lsm_h,                                                     &
+               surf_usm_h
 
 #if defined( __parallel )  &&  !defined( __mpifh )
@@ -1402,5 +1409,5 @@
 !--                            Determine surface level. Therefore, check for
 !--                            upward-facing wall on w-grid. 
-                               k_surf = get_topography_top_index_ji( jp, ip, 'w' )
+                               k_surf = topo_top_ind(jp,ip,3)
                                IF ( seed_follows_topography )  THEN
 !
@@ -3469,5 +3476,5 @@
 !--          above topography (Prandtl-layer height)
 !--          Determine vertical index of topography top
-             k_wall = get_topography_top_index_ji( jp, ip, 's' )
+             k_wall = topo_top_ind(jp,ip,0)
 
              IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
@@ -3558,5 +3565,5 @@
 !
 !--          Determine vertical index of topography top
-             k_wall = get_topography_top_index_ji( jp,ip, 's' )
+             k_wall = topo_top_ind(jp,ip,0)
 
              IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 4167)
+++ palm/trunk/SOURCE/modules.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +topo_top_ind
+! 
+! 4131 2019-08-02 11:06:18Z monakurppa
 ! Add max_pr_salsa to control_parameters. Used in creating profile output for
 ! salsa.
@@ -1755,4 +1758,5 @@
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  advc_flags_m            !< flags used to degrade order of advection scheme for momentum
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  advc_flags_s            !< flags used to degrade order of advection scheme for scalar quantities
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo_top_ind            !< precalculated topography top indices
     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  wall_flags_0            !< flags to mask topography and surface-bounded grid points
 
Index: palm/trunk/SOURCE/multi_agent_system_mod.f90
===================================================================
--- palm/trunk/SOURCE/multi_agent_system_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/multi_agent_system_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 3987 2019-05-22 09:52:13Z kanani
 ! Introduce alternative switch for debug output during timestepping
 ! 
@@ -131,4 +134,5 @@
     USE indices,                                                               &
         ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,        &
+               topo_top_ind,                                                   &
                wall_flags_0
 
@@ -374,6 +378,8 @@
  SUBROUTINE multi_agent_system
 
-    USE biometeorology_mod,                                                   &
-        ONLY:  bio_calc_ipt, bio_calculate_mrt_grid, bio_get_thermal_index_input_ij
+    USE biometeorology_mod,                                                    &
+        ONLY:  bio_calc_ipt,                                                   &
+               bio_calculate_mrt_grid,                                         &
+               bio_get_thermal_index_input_ij
 
 
@@ -3018,7 +3024,4 @@
            ONLY:  coupling_char, initializing_actions, io_blocks, io_group
 
-       USE surface_mod,                                                        &
-           ONLY:  get_topography_top_index, get_topography_top_index_ji
-
        USE arrays_3d,                                                          &
            ONLY:  zu, zw
@@ -3077,6 +3080,6 @@
        DO il = nxlg, nxrg
           DO jl = nysg, nyng
-             top_top_s(jl,il) = get_topography_top_index_ji(jl,il,'s') + 1
-             top_top_w(jl,il) = get_topography_top_index_ji(jl,il,'w')
+             top_top_s(jl,il) = topo_top_ind(jl,il,0) + 1
+             top_top_w(jl,il) = topo_top_ind(jl,il,3)
           ENDDO
        ENDDO
Index: palm/trunk/SOURCE/nesting_offl_mod.f90
===================================================================
--- palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/nesting_offl_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4125 2019-07-29 13:31:44Z suehring
 ! In order to enable netcdf parallel access, allocate dummy arrays for the 
 ! lateral boundary data on cores that actually do not belong to these 
@@ -119,5 +122,7 @@
     USE indices,                                                               &
         ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys,                  &
-               nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
+               nysv, nysg, nyn, nyng, nzb, nz, nzt,                            &
+               topo_top_ind,                                                   &
+               wall_flags_0
 
     USE kinds
@@ -873,7 +878,4 @@
        
        USE kinds
-       
-       USE surface_mod,                                                        &
-           ONLY:  get_topography_top_index, get_topography_top_index_ji
 
        IMPLICIT NONE 
@@ -918,5 +920,5 @@
 !
 !--          Determine topography top index at current (j,i) index
-             k_surface = get_topography_top_index_ji( j, i, 's' )
+             k_surface = topo_top_ind(j,i,0)
 !
 !--          Pre-compute surface virtual temperature. Therefore, use 2nd 
@@ -978,5 +980,5 @@
        
           DO  i = nxl, nxr
-             k_surface = get_topography_top_index_ji( j, i, 's' )
+             k_surface = topo_top_ind(j,i,0)
  
              IF ( humidity )  THEN
@@ -1032,5 +1034,5 @@
 !--    turbulence generator accordingly. If Rayleigh damping would be applied
 !--    near buildings, etc., this would spoil the simulation results. 
-       topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
+       topo_max_l = zw(MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ))
        
 #if defined( __parallel )
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4168)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4127 2019-07-30 14:47:10Z suehring
 ! Output of 3D plant canopy variables changed. It is now relative to the local
 ! terrain rather than located at the acutal vertical level in the model. This 
@@ -269,12 +272,9 @@
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
-               nz, nzb, nzt
+               nz, nzb, nzt, topo_top_ind
 
     USE kinds
 
     USE pegrid
-
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
 
 
@@ -1292,5 +1292,5 @@
 !--          Check whether topography and local vegetation on top exceed 
 !--          height of the model domain.
-             k = get_topography_top_index_ji( j, i, 's' )
+             k = topo_top_ind(j,i,0)
              IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
                 message_string =  'Local vegetation height on top of ' //      &
@@ -1591,5 +1591,5 @@
                    CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
                 ENDIF
-                kk = get_topography_top_index_ji( j, i, 's' )
+                kk = topo_top_ind(j,i,0)
                 lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
              CASE DEFAULT
@@ -1677,5 +1677,5 @@
 !
 !--                Determine topography-top index on u-grid
-                   k_wall = get_topography_top_index_ji( j, i, 'u' )
+                   k_wall = topo_top_ind(j,i,1)
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
 
@@ -1741,5 +1741,5 @@
 !
 !--                Determine topography-top index on v-grid
-                   k_wall = get_topography_top_index_ji( j, i, 'v' )
+                   k_wall = topo_top_ind(j,i,2)
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1806,5 +1806,5 @@
 !
 !--                Determine topography-top index on w-grid
-                   k_wall = get_topography_top_index_ji( j, i, 'w' )
+                   k_wall = topo_top_ind(j,i,3)
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
@@ -1858,5 +1858,5 @@
                    DO  j = nys, nyn
 !--                   Determine topography-top index on scalar-grid
-                      k_wall = get_topography_top_index_ji( j, i, 's' )
+                      k_wall = topo_top_ind(j,i,0)
                       DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
                          kk = k - k_wall   !- lad arrays are defined flat
@@ -1869,5 +1869,5 @@
                    DO  j = nys, nyn
 !--                   Determine topography-top index on scalar-grid
-                      k_wall = get_topography_top_index_ji( j, i, 's' )
+                      k_wall = topo_top_ind(j,i,0)
                       DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
                          kk = k - k_wall   !- lad arrays are defined flat
@@ -1885,5 +1885,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index_ji( j, i, 's' )
+                   k_wall = topo_top_ind(j,i,0)
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1923,5 +1923,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index_ji( j, i, 's' )
+                   k_wall = topo_top_ind(j,i,0)
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -1952,5 +1952,5 @@
 !
 !--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index_ji( j, i, 's' )
+                   k_wall = topo_top_ind(j,i,0)
 
                    DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
@@ -2046,5 +2046,5 @@
 !
 !--          Determine topography-top index on u-grid
-             k_wall = get_topography_top_index_ji( j, i, 'u' )
+             k_wall = topo_top_ind(j,i,1)
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
@@ -2106,5 +2106,5 @@
 !
 !--          Determine topography-top index on v-grid
-             k_wall = get_topography_top_index_ji( j, i, 'v' )
+             k_wall = topo_top_ind(j,i,2)
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -2166,5 +2166,5 @@
 !
 !--          Determine topography-top index on w-grid
-             k_wall = get_topography_top_index_ji( j, i, 'w' )
+             k_wall = topo_top_ind(j,i,3)
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
@@ -2213,5 +2213,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index_ji( j, i, 's' )
+             k_wall = topo_top_ind(j,i,0)
 
              IF ( humidity ) THEN
@@ -2233,5 +2233,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index_ji( j, i, 's' )
+             k_wall = topo_top_ind(j,i,0)
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -2266,5 +2266,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index_ji( j, i, 's' )
+             k_wall = topo_top_ind(j,i,0)
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
@@ -2292,5 +2292,5 @@
 !
 !--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index_ji( j, i, 's' )
+             k_wall = topo_top_ind(j,i,0)
 
              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 4167)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4029 2019-06-14 14:04:35Z raasch
 ! nest_chemistry switch removed
 ! 
@@ -478,5 +481,5 @@
     USE indices,                                                               &
         ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
-               nysv, nz, nzb, nzt, wall_flags_0
+               nysv, nz, nzb, nzt, topo_top_ind, wall_flags_0
 
     USE bulk_cloud_model_mod,                                                  &
@@ -528,5 +531,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE
@@ -2379,5 +2382,5 @@
           DO  j = nys, nyn
              sub_sum = 0.0_wp
-             k_wall = get_topography_top_index_ji( j, i, 'u' )
+             k_wall = topo_top_ind(j,i,1)
              DO   k = k_wall + 1, nzt
                 sub_sum = sub_sum + dzw(k)
@@ -2401,5 +2404,5 @@
           DO  j = nys, nyn
              sub_sum = 0.0_wp
-             k_wall = get_topography_top_index_ji( j, i, 'u' )
+             k_wall = topo_top_ind(j,i,1)
              DO   k = k_wall + 1, nzt
                 sub_sum = sub_sum + dzw(k)
@@ -2423,5 +2426,5 @@
           DO  i = nxl, nxr
              sub_sum = 0.0_wp
-             k_wall = get_topography_top_index_ji( j, i, 'v' )
+             k_wall = topo_top_ind(j,i,2)
              DO  k = k_wall + 1, nzt
                 sub_sum = sub_sum + dzw(k)
@@ -2445,5 +2448,5 @@
           DO  i = nxl, nxr
              sub_sum = 0.0_wp
-             k_wall = get_topography_top_index_ji( j, i, 'v' )
+             k_wall = topo_top_ind(j,i,2)
              DO  k = k_wall + 1, nzt
                 sub_sum = sub_sum + dzw(k)
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4168)
@@ -28,6 +28,5 @@
 ! -----------------
 ! $Id$
-! Changed behaviour of masked output over surface to follow terrain and ignore 
-! buildings (J.Resler, T.Gronemeier)
+! Replace function get_topography_top_index by topo_top_ind
 ! 
 ! 4157 2019-08-14 09:19:12Z suehring
@@ -691,5 +690,5 @@
     USE indices,                                                               &
         ONLY:  nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,   &
-               nzb, nzt
+               nzb, nzt, topo_top_ind
 
     USE, INTRINSIC :: iso_c_binding
@@ -737,6 +736,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
-               ind_pav_green, ind_veg_wall, ind_wat_win,                       &
+        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
                surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
                vertical_surfaces_exist
@@ -745,6 +743,4 @@
 
     CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
-
-    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
 
 !
@@ -3583,5 +3579,5 @@
 !
 !--       Determine minimum topography top index. 
-          k_topo_l = MINVAL( get_topography_top_index( 's' ) )
+          k_topo_l = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 #if defined( __parallel )
           CALL MPI_ALLREDUCE( k_topo_l, k_topo, 1, MPI_INTEGER, MPI_MIN, &
@@ -3773,5 +3769,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   k_topo_l = get_topography_top_index_ji( j, i, 's' )
+                   k_topo_l = topo_top_ind(j,i,0)
                    DO k = k_topo_l+1, nzt+1
                       rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k-k_topo_l)  * d_hours_day
@@ -4066,5 +4062,5 @@
 !
 !--             Obtain topography top index (lower bound of RRTMG)
-                k_topo = get_topography_top_index_ji( j, i, 's' )
+                k_topo = topo_top_ind(j,i,0)
 
                 IF ( lw_radiation )  THEN
@@ -5688,5 +5684,5 @@
 !
 !--          Following expression equals former kk = k - nzb_s_inner(j,i)
-             kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
+             kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
              pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
                  * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
@@ -5701,5 +5697,5 @@
                  j = pcbl(iy, ipcgb)
                  k = pcbl(iz, ipcgb)
-                 kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
+                 kk = k - topo_top_ind(j,i,0)  !- lad arrays are defined flat
                  CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), &
                                                    pc_transpiration_rate(kk,j,i), pc_latent_rate(kk,j,i) )
@@ -6217,6 +6213,6 @@
 !--    removed later). The following contruct finds the lowest / largest index
 !--    for any upward-facing wall (see bit 12).
-       nzubl = MINVAL( get_topography_top_index( 's' ) )
-       nzutl = MAXVAL( get_topography_top_index( 's' ) )
+       nzubl = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
+       nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
 
        nzubl = MAX( nzubl, nzb )
@@ -6235,5 +6231,5 @@
 !
 !--                Find topography top index
-                   k_topo = get_topography_top_index_ji( j, i, 's' )
+                   k_topo = topo_top_ind(j,i,0)
 
                    DO k = nzt+1, 0, -1
@@ -6358,5 +6354,5 @@
 !
 !--                Find topography top index
-                   k_topo = get_topography_top_index_ji( j, i, 's' )
+                   k_topo = topo_top_ind(j,i,0)
 
                    DO k = k_topo + 1, pct(j,i)
@@ -6793,5 +6789,5 @@
         ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
         nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
-        nzterrl = get_topography_top_index( 's' )
+        nzterrl = topo_top_ind(nys:nyn,nxl:nxr,0)
         CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
                             nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
@@ -6803,5 +6799,5 @@
         DEALLOCATE(nzterrl_l)
 #else
-        nzterr = RESHAPE( get_topography_top_index( 's' ), (/(nx+1)*(ny+1)/) )
+        nzterr = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
 #endif
         IF ( plant_canopy )  THEN
@@ -6884,5 +6880,5 @@
             DO i = nxl, nxr
                 DO j = nys, nyn
-                    k = get_topography_top_index_ji( j, i, 's' )
+                    k = topo_top_ind(j,i,0)
 
                     sub_lad(k:nz_plant_t, j, i) = lad_s(0:nz_plant_t-k, j, i)
@@ -10143,4 +10139,6 @@
     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
 
+    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
+
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
@@ -10540,4 +10538,6 @@
     LOGICAL      ::  found       !<
 
+    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
+
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
@@ -11107,13 +11107,12 @@
     CHARACTER (LEN=*) ::  variable   !< 
 
+    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
+
     INTEGER(iwp) ::  av              !< 
     INTEGER(iwp) ::  i               !< 
     INTEGER(iwp) ::  j               !< 
-    INTEGER(iwp) ::  k               !<
-    INTEGER(iwp) ::  im              !< loop index for masked variables
-    INTEGER(iwp) ::  jm              !< loop index for masked variables
-    INTEGER(iwp) ::  kk              !< 
+    INTEGER(iwp) ::  k               !< 
     INTEGER(iwp) ::  mid             !< masked output running index
-    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
+    INTEGER(iwp) ::  topo_top_index  !< k index of highest horizontal surface
 
     LOGICAL ::  found                !< true if output array was found
@@ -11127,6 +11126,8 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
 
+
+    found    = .TRUE.
+    grid     = 's'
     resorted = .FALSE.
-    found = .TRUE.
 
     SELECT CASE ( TRIM( variable ) )
@@ -11215,18 +11216,16 @@
              DO  j = 1, mask_size_l(mid,2)
 !
-!--             Get k index of the highest terraing surface
-                im = mask_i(mid,i)
-                jm = mask_j(mid,j)
-                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), &
-                              DIM = 1 ) - 1
+!--             Get k index of highest horizontal surface
+                topo_top_index = topo_top_ind(mask_j(mid,j), &
+                                              mask_i(mid,i),   &
+                                              0 )
+!
+!--             Save output array
                 DO  k = 1, mask_size_l(mid,3)
-                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
-!
-!--                Set value if not in building
-                   IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
-                      local_pf(i,j,k) = fill_value
-                   ELSE
-                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
-                   ENDIF
+                   local_pf(i,j,k) = to_be_resorted(                         &
+                                          MIN( topo_top_index+mask_k(mid,k), &
+                                               nzt+1 ),                      &
+                                          mask_j(mid,j),                     &
+                                          mask_i(mid,i)                     )
                 ENDDO
              ENDDO
Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 4168)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Remove functions get_topography_top_index. These are now replaced by 
+! precalculated arrays because of too much CPU-time consumption
+! 
+! 4159 2019-08-15 13:31:35Z suehring
 ! Surface classification revised and adjusted to changes in init_grid
 ! 
@@ -646,5 +650,5 @@
 
     PRIVATE
-
+    
     INTERFACE init_bc
        MODULE PROCEDURE init_bc
@@ -698,7 +702,5 @@
 !
 !-- Public subroutines and functions
-    PUBLIC get_topography_top_index,                                                               &
-           get_topography_top_index_ji,                                                            &
-           init_bc,                                                                                &
+    PUBLIC init_bc,                                                                                &
            init_single_surface_properties,                                                         &
            init_surfaces,                                                                          &
@@ -710,5 +712,6 @@
 
 #if defined( _OPENACC )
-    PUBLIC enter_surface_arrays, exit_surface_arrays
+    PUBLIC enter_surface_arrays,                                                                   &
+           exit_surface_arrays
 #endif
 
@@ -3141,90 +3144,4 @@
        
     END SUBROUTINE init_single_surface_properties
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Determines topography-top index at given (j,i)-position.  
-!------------------------------------------------------------------------------!
-    FUNCTION get_topography_top_index_ji( j, i, grid )
-
-       IMPLICIT NONE 
-
-       CHARACTER(LEN=*) ::  grid                         !< flag to distinquish between staggered grids
-       INTEGER(iwp)     ::  i                            !< grid index in x-dimension
-       INTEGER(iwp)     ::  ibit                         !< bit position where topography information is stored on respective grid
-       INTEGER(iwp)     ::  j                            !< grid index in y-dimension
-       INTEGER(iwp)     ::  get_topography_top_index_ji  !< topography top index
-
-       SELECT CASE ( TRIM( grid ) )
-
-          CASE ( 's'     )
-             ibit = 12
-          CASE ( 'u'     )
-             ibit = 14
-          CASE ( 'v'     )
-             ibit = 16
-          CASE ( 'w'     )
-             ibit = 18
-          CASE ( 's_out' )
-             ibit = 24
-          CASE DEFAULT
-!
-!--          Set default to scalar grid
-             ibit = 12 
-
-       END SELECT
-
-       get_topography_top_index_ji = MAXLOC(                                   &
-                                     MERGE( 1, 0,                              &
-                                            BTEST( wall_flags_0(:,j,i), ibit ) &
-                                          ), DIM = 1                           &
-                                           ) - 1
-
-       RETURN
-
-    END FUNCTION get_topography_top_index_ji
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Determines topography-top index at each (j,i)-position.  
-!------------------------------------------------------------------------------!
-    FUNCTION get_topography_top_index( grid )
-
-       IMPLICIT NONE 
-
-       CHARACTER(LEN=*) ::  grid                      !< flag to distinquish between staggered grids
-       INTEGER(iwp)     ::  ibit                      !< bit position where topography information is stored on respective grid
-       INTEGER(iwp), DIMENSION(nys:nyn,nxl:nxr) ::  get_topography_top_index  !< topography top index
-
-       SELECT CASE ( TRIM( grid ) )
-
-          CASE ( 's'     )
-             ibit = 12
-          CASE ( 'u'     )
-             ibit = 14
-          CASE ( 'v'     )
-             ibit = 16
-          CASE ( 'w'     )
-             ibit = 18
-          CASE ( 's_out' )
-             ibit = 24
-          CASE DEFAULT
-!
-!--          Set default to scalar grid
-             ibit = 12 
-
-       END SELECT
-
-       get_topography_top_index(nys:nyn,nxl:nxr) = MAXLOC(                     &
-                         MERGE( 1, 0,                                          &
-                                 BTEST( wall_flags_0(:,nys:nyn,nxl:nxr), ibit )&
-                              ), DIM = 1                                       &
-                                                         ) - 1
-
-       RETURN
-
-    END FUNCTION get_topography_top_index
 
 !------------------------------------------------------------------------------!
Index: palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4110 2019-07-22 17:05:21Z suehring
 ! pass integer flag array as well as boundary flags to WS scalar advection 
 ! routine
@@ -254,4 +257,5 @@
         ONLY:  advc_flags_s,                                                  &
                nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,    &
+               topo_top_ind,                                                  &
                wall_flags_0
 
@@ -272,5 +276,4 @@
         ONLY:  bc_h,                                                          &
                bc_v,                                                          &
-               get_topography_top_index_ji,                                   &
                surf_def_h,                                                    &
                surf_def_v,                                                    &
@@ -1365,5 +1368,5 @@
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
+                nz_s_shift = topo_top_ind(j,i,0)
 
                 e(nz_s_shift:nzt+1,j,i)  =  e(0:nzt+1-nz_s_shift,j,i)
@@ -1375,5 +1378,5 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   nz_s_shift = get_topography_top_index_ji( j, i, 's' )
+                   nz_s_shift = topo_top_ind(j,i,0)
 
                    diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
@@ -1394,5 +1397,5 @@
              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
                   nysg <= 0 .AND. nyng >= 0        )  THEN
-                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
+                nz_s_shift_l = topo_top_ind(0,0,0)
              ELSE
                 nz_s_shift_l = 0
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4168)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 4148 2019-08-08 11:26:00Z suehring
 ! - Add anthropogenic heat output factors for heating and cooling to building
 !   data base
@@ -538,5 +541,5 @@
     USE indices,                                                               &
         ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
-               nysg, nzb, nzt, nbgp, wall_flags_0
+               nysg, nzb, nzt, nbgp, topo_top_ind, wall_flags_0
 
     USE, INTRINSIC :: iso_c_binding 
@@ -559,6 +562,5 @@
 
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji, get_topography_top_index,          &
-               ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h,           &
+        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h,           &
                surf_usm_v, surface_restore_elements
 
@@ -6159,5 +6161,5 @@
                     READ( 151, *, err=12, end=13 )  i, j, k, heat
                     IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
-                        IF ( k <= naheatlayers  .AND.  k > get_topography_top_index_ji( j, i, 's' ) )  THEN
+                        IF ( k <= naheatlayers  .AND.  k > topo_top_ind(j,i,0) )  THEN
 !--                         write heat into the array
                             aheat(k,j,i) = heat
@@ -7394,6 +7396,6 @@
               ELSE
                  WRITE(9,*) 'Problem reading USM data:'
-                 WRITE(9,*) l,i,j,kw,get_topography_top_index_ji( j, i, 's' )
-                 WRITE(9,*) ii,iw,jw,kw,get_topography_top_index_ji( jw, iw, 's' )
+                 WRITE(9,*) l,i,j,kw,topo_top_ind(j,i,0)
+                 WRITE(9,*) ii,iw,jw,kw,topo_top_ind(jw,iw,0)
                  WRITE(9,*) usm_par(ii,jw,iw),usm_par(ii+1,jw,iw)
                  WRITE(9,*) usm_par(ii+2,jw,iw),usm_par(ii+3,jw,iw)
@@ -8614,9 +8616,9 @@
 !--         TO_DO: activate, if testcase is available
 !--         !$OMP PARALLEL DO PRIVATE (i, j, k, acoef, rho_cp)
-!--         it may also improve performance to move get_topography_top_index_ji before the k-loop
+!--         it may also improve performance to move topo_top_ind before the k-loop
             DO i = nxl, nxr
                DO j = nys, nyn
                   DO k = nz_urban_b, min(nz_urban_t,naheatlayers)
-                     IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN
+                     IF ( k > topo_top_ind(j,i,0) ) THEN
 !
 !--                    increase of pt in box i,j,k in time dt_3d 
Index: palm/trunk/SOURCE/user_data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/user_data_output_mask.f90	(revision 4167)
+++ palm/trunk/SOURCE/user_data_output_mask.f90	(revision 4168)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Remove dependency on surface_mod + example for terrain-following output 
+! adjusted
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration
@@ -72,7 +76,4 @@
     USE kinds
     
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
-    
     USE user
 
@@ -80,5 +81,4 @@
 
     CHARACTER (LEN=*) ::  variable  !< 
-    CHARACTER (LEN=5) ::  grid      !< flag to distinquish between staggered grids
 
     INTEGER(iwp) ::  av             !< 
@@ -87,5 +87,5 @@
 !    INTEGER(iwp) ::  j              !<
 !    INTEGER(iwp) ::  k              !<
-!    INTEGER(iwp) ::  topo_top_ind   !< k index of highest horizontal surface
+!    INTEGER(iwp) ::  topo_top_index !< k index of highest horizontal surface
 
     LOGICAL ::  found               !< 
@@ -102,5 +102,4 @@
 
     found = .TRUE.
-    grid  = 's'
 
     SELECT CASE ( TRIM( variable ) )
@@ -131,14 +130,14 @@
 !!
 !!--                   Get k index of highest horizontal surface
-!                      topo_top_ind = get_topography_top_index_ji( &
+!                      topo_top_index = topo_top_ind( &
 !                                        mask_j(mid,j), &
 !                                        mask_i(mid,i), &
-!                                        grid )
+!                                        1          )
 !!
 !!--                   Save output array
 !                      DO  k = 1, mask_size_l(mid,3)
-!                         local_pf(i,j,k) = u2(MIN( topo_top_ind+mask_k(mid,k),&
-!                                                   nzt+1 ),                   &
-!                                              mask_j(mid,j),                  &
+!                         local_pf(i,j,k) = u2(MIN( topo_top_index+mask_k(mid,k),&
+!                                                   nzt+1 ),                     &
+!                                              mask_j(mid,j),                    &
 !                                              mask_i(mid,i)                   )
 !                      ENDDO
@@ -166,15 +165,15 @@
 !!
 !!--                   Get k index of highest horizontal surface
-!                      topo_top_ind = get_topography_top_index_ji( &
+!                      topo_top_index = topo_top_ind(   &
 !                                        mask_j(mid,j), &
 !                                        mask_i(mid,i), &
-!                                        grid )
+!                                        1 )
 !!
 !!--                   Save output array
 !                      DO  k = 1, mask_size_l(mid,3)
-!                         local_pf(i,j,k) = u2_av(                             &
-!                                              MIN( topo_top_ind+mask_k(mid,k),&
-!                                                   nzt+1 ),                   &
-!                                              mask_j(mid,j),                  &
+!                         local_pf(i,j,k) = u2_av(                               &
+!                                              MIN( topo_top_index+mask_k(mid,k),&
+!                                                   nzt+1 ),                     &
+!                                              mask_j(mid,j),                    &
 !                                              mask_i(mid,i)                   )
 !                      ENDDO
Index: palm/trunk/SOURCE/virtual_measurement_mod.f90
===================================================================
--- palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4167)
+++ palm/trunk/SOURCE/virtual_measurement_mod.f90	(revision 4168)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Replace function get_topography_top_index by topo_top_ind
+! 
+! 3988 2019-05-22 11:32:37Z kanani
 ! Add variables to enable steering of output interval for virtual measurements
 ! 
@@ -115,5 +118,5 @@
 
     USE indices,                                                               &
-        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, wall_flags_0
+        ONLY:  nzb, nzt, nxl, nxr, nys, nyn, nx, ny, topo_top_ind, wall_flags_0
 
     USE kinds
@@ -435,7 +438,4 @@
                netcdf_data_input_get_dimension_length,                         &
                netcdf_data_input_att, netcdf_data_input_var
-               
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
         
     IMPLICIT NONE
@@ -763,5 +763,5 @@
 !--                Determine vertical index which correspond to the observation
 !--                height.
-                   ksurf = get_topography_top_index_ji( js, is, 's' )
+                   ksurf = topo_top_ind(js,is,0)
                    ks = MINLOC( ABS( zu - zw(ksurf) - z_ag(t,n) ), DIM = 1 ) - 1
 !
@@ -822,6 +822,5 @@
                       vmea(l)%j(ns) = j
                       vmea(l)%k(ns) = k
-                      vmea(l)%z_ag(ns)  = zu(k) -                              &
-                                   zw(get_topography_top_index_ji( j, i, 's' ))
+                      vmea(l)%z_ag(ns)  = zu(k) - zw(topo_top_ind(j,i,0))
                    ENDIF
                 ENDDO
