Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 2967)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 2968)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix in output of timeseries quantities in case of elevated model surfaces.
+! 
+! 2817 2018-02-19 16:32:21Z knoop
 ! Preliminary gust module interface implemented
 ! 
@@ -292,5 +295,6 @@
     USE indices,                                                               &
         ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
-                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
+                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, topo_min_level,     &
+                wall_flags_0
         
     USE kinds
@@ -2102,26 +2106,30 @@
 
 !
-!--    Collect the time series quantities
-       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
-       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
+!--    Collect the time series quantities. Please note, timeseries quantities 
+!--    which are collected from horizontally averaged profiles, e.g. wpt 
+!--    or pt(zp), are treated specially. In case of elevated model surfaces, 
+!--    index nzb+1 might be within topography and data will be zero. Therefore,
+!--    take value for the first atmosphere index, which is topo_min_level+1. 
+       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
+       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
        ts_value(3,sr) = dt_3d
-       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
-       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
+       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)          ! u*
+       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)        ! th*
        ts_value(6,sr) = u_max
        ts_value(7,sr) = v_max
        ts_value(8,sr) = w_max
-       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
-       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
-       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
-       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
-       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
-       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
-       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
-       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
-       ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)   ! pt(0)
-       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
-       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
-       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
-       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
+       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)       ! new divergence
+       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)       ! old Divergence
+       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)       ! z_i(1)
+       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)       ! z_i(2)
+       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)       ! w*
+       ts_value(14,sr) = hom(nzb,1,16,sr)              ! w'pt'   at k=0
+       ts_value(15,sr) = hom(topo_min_level+1,1,16,sr) ! w'pt'   at k=1
+       ts_value(16,sr) = hom(topo_min_level+1,1,18,sr) ! wpt     at k=1
+       ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)      ! pt(0)
+       ts_value(18,sr) = hom(topo_min_level+1,1,4,sr)  ! pt(zp)
+       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)       ! u'w'    at k=0
+       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)       ! v'w'    at k=0
+       ts_value(21,sr) = hom(nzb,1,48,sr)              ! w"q"    at k=0
 
        IF ( .NOT. neutral )  THEN
Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 2967)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 2968)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfix in initialization in case of elevated model surface. Introduce
+! index for minimum topography-top. 
+! 
+! 2955 2018-04-09 15:14:01Z suehring
 ! Improve topography filter routine and add ghost-point exchange for building 
 ! ID and building type. 
@@ -332,5 +336,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
+               nzb_w_outer, nzt, topo_min_level
     
     USE kinds
@@ -348,11 +352,11 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i             !< index variable along x 
-    INTEGER(iwp) ::  j             !< index variable along y
-    INTEGER(iwp) ::  k             !< index variable along z
-    INTEGER(iwp) ::  k_top         !< topography top index on local PE
-    INTEGER(iwp) ::  l             !< loop variable 
-    INTEGER(iwp) ::  nzb_local_max !< vertical grid index of maximum topography height
-    INTEGER(iwp) ::  nzb_local_min !< vertical grid index of minimum topography height
+    INTEGER(iwp) ::  i                !< index variable along x 
+    INTEGER(iwp) ::  j                !< index variable along y
+    INTEGER(iwp) ::  k                !< index variable along z
+    INTEGER(iwp) ::  k_top            !< topography top index on local PE
+    INTEGER(iwp) ::  l                !< loop variable 
+    INTEGER(iwp) ::  nzb_local_max    !< vertical grid index of maximum topography height
+    INTEGER(iwp) ::  nzb_local_min    !< vertical grid index of minimum topography height
                                       
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
@@ -553,10 +557,20 @@
 !   
 !-- Finally, if topography extents up to the model top, limit nzb_max to nzt.
-    nzb_max = MIN( nzb_max, nzt ) 
-
+    nzb_max = MIN( nzb_max, nzt )
+!
+!-- Determine minimum index of topography. Usually, this will be nzb. In case
+!-- there is elevated topography, however, the lowest topography will be higher. 
+!-- This index is e.g. used to calculate mean first-grid point atmosphere 
+!-- temperature, surface pressure and density, etc. .
+    topo_min_level   = 0
+#if defined( __parallel )
+    CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),             &
+                        topo_min_level, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
+#else
+    topo_min_level = MINVAL( get_topography_top_index( 's' ) )
+#endif
 !
 !-- Initialize boundary conditions via surface type
     CALL init_bc
-
 !
 !-- Allocate and set topography height arrays required for data output
@@ -621,15 +635,16 @@
     nzb_local(nys:nyn,nxl:nxr) = get_topography_top_index( 's' )
     CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
-
+!
+!-- Check topography for consistency with model domain. Therefore, use
+!-- maximum and minium topography-top indices. Note, minimum topography top
+!-- index is already calculated.  
     IF ( TRIM( topography ) /= 'flat' )  THEN
 #if defined( __parallel )
        CALL MPI_ALLREDUCE( MAXVAL( get_topography_top_index( 's' ) ),          &
-                           nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
-       CALL MPI_ALLREDUCE( MINVAL( get_topography_top_index( 's' ) ),          &
-                           nzb_local_min, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )                   
+                           nzb_local_max, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )             
 #else
        nzb_local_max = MAXVAL( get_topography_top_index( 's' ) )
-       nzb_local_min = MINVAL( get_topography_top_index( 's' ) )
 #endif
+       nzb_local_min = topo_min_level
 !
 !--    Consistency checks
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2967)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2968)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix in initialization in case of elevated model surface
+! 
+! 2963 2018-04-12 14:47:44Z suehring
 ! - In initialization of surface types, consider the case that surface_fractions
 !   is not given in static input file.
@@ -2374,5 +2377,5 @@
 
        USE indices,                                                            &
-           ONLY:  nx, ny
+           ONLY:  nx, ny, topo_min_level
 
        USE pmc_interface,                                                      &
@@ -2408,5 +2411,5 @@
        IF (  .NOT.  cloud_physics )  THEN
           CALL calc_mean_profile( pt, 4 )
-          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn )
+          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(topo_min_level+1,1,4,0) * exn )
        ENDIF
 
@@ -2424,6 +2427,6 @@
        tt_soil_h_m%var_2d    = 0.0_wp
        tm_soil_h_m%var_2d    = 0.0_wp
-       tm_liq_h_m%var_1d  = 0.0_wp
-       surf_lsm_h%c_liq = 0.0_wp
+       tm_liq_h_m%var_1d     = 0.0_wp
+       surf_lsm_h%c_liq      = 0.0_wp
 
        surf_lsm_h%ghf = 0.0_wp
@@ -2443,6 +2446,6 @@
           tt_soil_v_m(l)%var_2d    = 0.0_wp
           tm_soil_v_m(l)%var_2d    = 0.0_wp
-          tm_liq_v_m(l)%var_1d  = 0.0_wp
-          surf_lsm_v(l)%c_liq = 0.0_wp
+          tm_liq_v_m(l)%var_1d     = 0.0_wp
+          surf_lsm_v(l)%c_liq      = 0.0_wp
 
           surf_lsm_v(l)%ghf = 0.0_wp
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 2967)
+++ /palm/trunk/SOURCE/modules.f90	(revision 2968)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +topo_min_level
+! 
+! 2964 2018-04-12 16:04:03Z raasch
 ! *_time_count variables are all initialized with zero now
 ! 
@@ -1719,33 +1722,34 @@
     USE kinds
 
-    INTEGER(iwp) ::  nbgp = 3     !< number of boundary ghost points
-    INTEGER(iwp) ::  ngp_sums     !< number of vertical profile grid points time number of output profiles - used for allreduce statements in MPI calls
-    INTEGER(iwp) ::  ngp_sums_ls  !< number of vertical profile grid points time number of large-scale forcing profiles - used for allreduce statements in MPI calls
-    INTEGER(iwp) ::  nnx          !< number of subdomain grid points in x-direction
-    INTEGER(iwp) ::  nx = 0       !< nx+1 = total number of grid points in x-direction
-    INTEGER(iwp) ::  nx_a         !< in coupled atmosphere-ocean runs: total number of grid points along x (atmosphere)
-    INTEGER(iwp) ::  nx_o         !< in coupled atmosphere-ocean runs: total number of grid points along x (ocean)
-    INTEGER(iwp) ::  nxl          !< left-most grid index of subdomain (excluding ghost points)
-    INTEGER(iwp) ::  nxlg         !< left-most grid index of subdomain (including ghost points)
-    INTEGER(iwp) ::  nxlu         !< =nxl+1 (at left domain boundary with inflow from left), else =nxl (used for u-velocity component)
-    INTEGER(iwp) ::  nxr          !< right-most grid index of subdomain (excluding ghost points)
-    INTEGER(iwp) ::  nxrg         !< right-most grid index of subdomain (including ghost points)
-    INTEGER(iwp) ::  nx_on_file   !< nx of previous run in job chain
-    INTEGER(iwp) ::  nny          !< number of subdomain grid points in y-direction
-    INTEGER(iwp) ::  ny = 0       !< ny+1 = total number of grid points in y-direction
-    INTEGER(iwp) ::  ny_a         !< in coupled atmosphere-ocean runs: total number of grid points along y (atmosphere)
-    INTEGER(iwp) ::  ny_o         !< in coupled atmosphere-ocean runs: total number of grid points along y (ocean)
-    INTEGER(iwp) ::  nyn          !< north-most grid index of subdomain (excluding ghost points)
-    INTEGER(iwp) ::  nyng         !< north-most grid index of subdomain (including ghost points)
-    INTEGER(iwp) ::  nys          !< south-most grid index of subdomain (excluding ghost points)
-    INTEGER(iwp) ::  nysg         !< south-most grid index of subdomain (including ghost points)
-    INTEGER(iwp) ::  nysv         !< =nys+1 (at south domain boundary with inflow from south), else =nys (used for v-velocity component)
-    INTEGER(iwp) ::  ny_on_file   !< ny of previous run in job chain
-    INTEGER(iwp) ::  nnz          !< number of subdomain grid points in z-direction
-    INTEGER(iwp) ::  nz = 0       !< total number of grid points in z-direction
-    INTEGER(iwp) ::  nzb          !< bottom grid index of computational domain
-    INTEGER(iwp) ::  nzb_diff     !< will be removed
-    INTEGER(iwp) ::  nzb_max      !< vertical index of topography top
-    INTEGER(iwp) ::  nzt          !< nzt+1 = top grid index of computational domain
+    INTEGER(iwp) ::  nbgp = 3       !< number of boundary ghost points
+    INTEGER(iwp) ::  ngp_sums       !< number of vertical profile grid points time number of output profiles - used for allreduce statements in MPI calls
+    INTEGER(iwp) ::  ngp_sums_ls    !< number of vertical profile grid points time number of large-scale forcing profiles - used for allreduce statements in MPI calls
+    INTEGER(iwp) ::  nnx            !< number of subdomain grid points in x-direction
+    INTEGER(iwp) ::  nx = 0         !< nx+1 = total number of grid points in x-direction
+    INTEGER(iwp) ::  nx_a           !< in coupled atmosphere-ocean runs: total number of grid points along x (atmosphere)
+    INTEGER(iwp) ::  nx_o           !< in coupled atmosphere-ocean runs: total number of grid points along x (ocean)
+    INTEGER(iwp) ::  nxl            !< left-most grid index of subdomain (excluding ghost points)
+    INTEGER(iwp) ::  nxlg           !< left-most grid index of subdomain (including ghost points)
+    INTEGER(iwp) ::  nxlu           !< =nxl+1 (at left domain boundary with inflow from left), else =nxl (used for u-velocity component)
+    INTEGER(iwp) ::  nxr            !< right-most grid index of subdomain (excluding ghost points)
+    INTEGER(iwp) ::  nxrg           !< right-most grid index of subdomain (including ghost points)
+    INTEGER(iwp) ::  nx_on_file     !< nx of previous run in job chain
+    INTEGER(iwp) ::  nny            !< number of subdomain grid points in y-direction
+    INTEGER(iwp) ::  ny = 0         !< ny+1 = total number of grid points in y-direction
+    INTEGER(iwp) ::  ny_a           !< in coupled atmosphere-ocean runs: total number of grid points along y (atmosphere)
+    INTEGER(iwp) ::  ny_o           !< in coupled atmosphere-ocean runs: total number of grid points along y (ocean)
+    INTEGER(iwp) ::  nyn            !< north-most grid index of subdomain (excluding ghost points)
+    INTEGER(iwp) ::  nyng           !< north-most grid index of subdomain (including ghost points)
+    INTEGER(iwp) ::  nys            !< south-most grid index of subdomain (excluding ghost points)
+    INTEGER(iwp) ::  nysg           !< south-most grid index of subdomain (including ghost points)
+    INTEGER(iwp) ::  nysv           !< =nys+1 (at south domain boundary with inflow from south), else =nys (used for v-velocity component)
+    INTEGER(iwp) ::  ny_on_file     !< ny of previous run in job chain
+    INTEGER(iwp) ::  nnz            !< number of subdomain grid points in z-direction
+    INTEGER(iwp) ::  nz = 0         !< total number of grid points in z-direction
+    INTEGER(iwp) ::  nzb            !< bottom grid index of computational domain
+    INTEGER(iwp) ::  nzb_diff       !< will be removed
+    INTEGER(iwp) ::  nzb_max        !< vertical index of topography top
+    INTEGER(iwp) ::  nzt            !< nzt+1 = top grid index of computational domain
+    INTEGER(iwp) ::  topo_min_level !< minimum topography-top index (usually equal to nzb)
 
     INTEGER(idp), DIMENSION(:), ALLOCATABLE ::  ngp_3d        !< number of grid points of the total domain
