Index: /palm/trunk/SOURCE/calc_mean_profile.f90
===================================================================
--- /palm/trunk/SOURCE/calc_mean_profile.f90	(revision 1737)
+++ /palm/trunk/SOURCE/calc_mean_profile.f90	(revision 1738)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! bugfix: if a layer is completely filled with topography, no mean is calculated
 ! 
 ! Former revisions:
@@ -126,5 +126,9 @@
 #endif
 
-          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
+          DO  k = nzb, nzt+1
+             IF ( ngp_2dh_s_inner(k,0) /= 0 )  THEN
+                hom(k,1,pr,0) = sums(k,pr) / ngp_2dh_s_inner(k,0)
+             ENDIF
+          ENDDO
 
        ENDIF
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 1737)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 1738)
@@ -19,5 +19,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! bugfixes for calculations in statistical regions which do not contain grid
+! points in the lowest vertical levels, mean surface level height considered
+! in the calculation of the characteristic vertical velocity,
+! old upstream parts removed
 ! 
 ! Former revisions:
@@ -222,4 +225,5 @@
     INTEGER(iwp) ::  j                   !<
     INTEGER(iwp) ::  k                   !<
+    INTEGER(iwp) ::  k_surface_level     !<
     INTEGER(iwp) ::  nt                  !<
     INTEGER(iwp) ::  omp_get_thread_num  !<
@@ -1247,27 +1251,29 @@
 !--    Final values are obtained by division by the total number of grid points 
 !--    used for summation. After that store profiles.
+!--    Check, if statistical regions do contain at least one grid point at the
+!--    respective k-level, otherwise division by zero will lead to undefined
+!--    values, which may cause e.g. problems with NetCDF output
 !--    Profiles:
        DO  k = nzb, nzt+1
-          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
-          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
-          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
-          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
-          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
-          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
-          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
-          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
-          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
-          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
-          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
-          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
-          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
-          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
-          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
-          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
-          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
+          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
+          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
+          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
+          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
+          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
+          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
+          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
+          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
+          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
+             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
+             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
+             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
+             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
+             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
+             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
+             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
+             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
+          ENDIF
        ENDDO
 
-!--    Upstream-parts
-       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
 !--    u* and so on
 !--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
@@ -1422,7 +1428,4 @@
        hom(:,1,114,sr) = sums(:,114)            !: L
 
-       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
-                                       ! upstream-parts u_x, u_y, u_z, v_x,
-                                       ! v_y, usw. (in last but one profile)
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -1445,11 +1448,9 @@
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp)  THEN
+             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
-                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
+             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5_wp * height )  THEN
@@ -1463,11 +1464,9 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
+             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               & 
-                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
+             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5_wp * height )  THEN
@@ -1519,17 +1518,26 @@
 
 !
+!--    Determine vertical index which is nearest to the mean surface level
+!--    height of the respective statistic region
+       DO  k = nzb, nzt
+          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
+             k_surface_level = k
+             EXIT
+          ENDIF
+       ENDDO
+!
 !--    Computation of both the characteristic vertical velocity and
 !--    the characteristic convective boundary layer temperature.
-!--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
-           .AND.  z_i(1) /= 0.0_wp )  THEN
-          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
-                                       hom(nzb,1,18,sr) *                      &
-                                       ABS( z_i(1) ) )**0.333333333_wp
-!--       so far this only works if Prandtl layer is used
-          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
+!--    The inversion height entering into the equation is defined with respect
+!--    to the mean surface level height of the respective statistic region.
+!--    The horizontal average at surface level index + 1 is input for the
+!--    average temperature.
+       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
+       THEN
+          hom(nzb+8,1,pr_palm,sr) = &
+             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
           hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
-          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
        ENDIF
 
@@ -1689,4 +1697,5 @@
     INTEGER(iwp) ::  j                   !<
     INTEGER(iwp) ::  k                   !<
+    INTEGER(iwp) ::  k_surface_level     !<
     INTEGER(iwp) ::  nt                  !<
     INTEGER(iwp) ::  omp_get_thread_num  !<
@@ -3230,27 +3239,29 @@
 !--    Final values are obtained by division by the total number of grid points
 !--    used for summation. After that store profiles.
+!--    Check, if statistical regions do contain at least one grid point at the
+!--    respective k-level, otherwise division by zero will lead to undefined
+!--    values, which may cause e.g. problems with NetCDF output
 !--    Profiles:
        DO  k = nzb, nzt+1
-          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
-          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
-          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
-          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
-          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
-          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
-          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
-          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
-          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
-          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
-          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
-          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
-          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
-          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
-          sums(k,89:114)           = sums(k,89:114)     / ngp_2dh(sr)
-          sums(k,115:pr_palm-2)    = sums(k,115:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
-
+          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
+          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
+          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
+          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
+          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
+          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
+          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
+          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
+          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
+             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
+             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
+             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
+             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
+             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
+             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
+             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
+             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
+          ENDIF
        ENDDO
 
-!--    Upstream-parts
-       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
 !--    u* and so on
 !--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
@@ -3271,7 +3282,9 @@
        IF ( max_pr_user > 0 )  THEN
           DO  k = nzb, nzt+1
-             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
-                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
-                                    ngp_2dh_s_inner(k,sr)
+             IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
+                sums(k,pr_palm+1:pr_palm+max_pr_user) = &
+                                       sums(k,pr_palm+1:pr_palm+max_pr_user) / &
+                                       ngp_2dh_s_inner(k,sr)
+             ENDIF
           ENDDO
        ENDIF
@@ -3363,7 +3376,4 @@
        END IF
 
-       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
-                                       ! upstream-parts u_x, u_y, u_z, v_x,
-                                       ! v_y, usw. (in last but one profile)
        hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
@@ -3386,11 +3396,9 @@
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
+             IF (  first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
-                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
+             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5_wp * height )  THEN
@@ -3404,11 +3412,9 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
-                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
+             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. &
-                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
+             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5_wp * height )  THEN
@@ -3460,17 +3466,26 @@
 
 !
+!--    Determine vertical index which is nearest to the mean surface level
+!--    height of the respective statistic region
+       DO  k = nzb, nzt
+          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
+             k_surface_level = k
+             EXIT
+          ENDIF
+       ENDDO
+
+!
 !--    Computation of both the characteristic vertical velocity and
 !--    the characteristic convective boundary layer temperature.
-!--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
-           .AND.  z_i(1) /= 0.0_wp )  THEN
-          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
-                                       hom(nzb,1,18,sr) *                      &
-                                       ABS( z_i(1) ) )**0.333333333_wp
-!--       so far this only works if Prandtl layer is used
-          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
+!--    The inversion height entering into the equation is defined with respect
+!--    to the mean surface level height of the respective statistic region.
+!--    The horizontal average at surface level index + 1 is input for the
+!--    average temperature.
+       IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
+          hom(nzb+8,1,pr_palm,sr) = &
+             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
+             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
        ELSE
           hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
-          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
        ENDIF
 
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 1737)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 1738)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! calculate mean surface level height for each statistic region
 !
 ! Former revisions:
@@ -273,8 +273,8 @@
     
     USE statistics,                                                            &
-        ONLY:  hom, hom_sum, pr_palm, rmask, spectrum_x, spectrum_y,           &
-               statistic_regions, sums, sums_divnew_l, sums_divold_l, sums_l,  &
-               sums_l_l, sums_up_fraction_l, sums_wsts_bc_l, ts_value,         &
-               var_d, weight_pres, weight_substep 
+        ONLY:  hom, hom_sum, mean_surface_level_height, pr_palm, rmask,        &
+               spectrum_x, spectrum_y, statistic_regions, sums, sums_divnew_l, &
+               sums_divold_l, sums_l, sums_l_l, sums_up_fraction_l,            &
+               sums_wsts_bc_l, ts_value, var_d, weight_pres, weight_substep
  
     USE surface_layer_fluxes_mod,                                              &
@@ -299,4 +299,5 @@
     REAL(wp), DIMENSION(1:2) ::  volume_flow_initial_l  !<
 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  mean_surface_level_height_l    !<
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l    !<
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_tmp  !<
@@ -306,5 +307,7 @@
 !
 !-- Allocate arrays
-    ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
+    ALLOCATE( mean_surface_level_height(0:statistic_regions),               &
+              mean_surface_level_height_l(0:statistic_regions),             &
+              ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
               ngp_3d(0:statistic_regions),                                  &
               ngp_3d_inner(0:statistic_regions),                            &
@@ -1820,4 +1823,5 @@
 !
 !-- Compute total sum of active mask grid points
+!-- and the mean surface level height for each statistic region
 !-- ngp_2dh: number of grid points of a horizontal cross section through the
 !--          total domain
@@ -1834,4 +1838,7 @@
     ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
 
+    mean_surface_level_height   = 0.0_wp
+    mean_surface_level_height_l = 0.0_wp
+
     DO  sr = 0, statistic_regions
        DO  i = nxl, nxr
@@ -1841,4 +1848,6 @@
 !--             All xy-grid points
                 ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
+                mean_surface_level_height_l(sr) = mean_surface_level_height_l(sr) + &
+                                                  zw(nzb_s_inner(j,i))
 !
 !--             xy-grid points above topography
@@ -1873,4 +1882,9 @@
                         MPI_SUM, comm2d, ierr )
     ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( mean_surface_level_height_l(0),                       &
+                        mean_surface_level_height(0), sr, MPI_REAL,           &
+                        MPI_SUM, comm2d, ierr )
+    mean_surface_level_height = mean_surface_level_height / REAL( ngp_2dh )
 #else
     ngp_2dh         = ngp_2dh_l
@@ -1878,4 +1892,5 @@
     ngp_2dh_s_inner = ngp_2dh_s_inner_l
     ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
+    mean_surface_level_height = mean_surface_level_height_l / REAL( ngp_2dh_l )
 #endif
 
@@ -1892,5 +1907,6 @@
     ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
 
-    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp )
+    DEALLOCATE( mean_surface_level_height_l, ngp_2dh_l, ngp_2dh_outer_l,      &
+                ngp_3d_inner_l, ngp_3d_inner_tmp )
 
     CALL location_message( 'leaving init_3d_model', .TRUE. )
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1737)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1738)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +mean_surface_level_height
 ! 
 ! Former revisions:
@@ -1467,5 +1467,6 @@
     LOGICAL ::  flow_statistics_called = .FALSE.
     REAL(wp) ::     u_max, v_max, w_max
-    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  sums_divnew_l, sums_divold_l,   &
+    REAL(wp), DIMENSION(:), ALLOCATABLE       ::  mean_surface_level_height,      &
+                                                  sums_divnew_l, sums_divold_l,   &
                                                   var_d, weight_substep,          &
                                                   weight_pres
