Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 4629)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 4630)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! In case of ASCII topography input flag grid points as terrain and building.
+!
+! 4601 2020-07-14 12:06:09Z suehring
 ! Minor formatting adjustments
 ! 
@@ -1347,8 +1350,11 @@
 !--             Flag topography for all grid points which are below
 !--             the local topography height. 
-!--             Note, each topography is flagged as building. 
+!--             Note, each topography is flagged as building (bit 2) as well as 
+!--             terrain (bit 1) in order to employ urban-surface as well as 
+!--             land-surface model.
                 IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) )  THEN
                    topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
-                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates building
+                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
+                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
                 ENDIF
              ENDDO
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4629)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 4630)
@@ -25,4 +25,11 @@
 ! -----------------
 ! $Id$
+! - Bugfix in level 3 initialization of root-area-density
+! - Bugfix in resistance calculation - avoid potential divisions by zero
+! - Avoid double classifiation of vertical walls (at surfaces that are also
+!   covered by buildings)
+! - Minor formatting adjustment to increase readability
+!
+! 4602 2020-07-14 14:49:45Z suehring
 ! - Bugfix in level 3 initialization of pavements - wrongly assumed existence of
 !   pavement_subsurface_pars
@@ -1579,6 +1586,4 @@
     ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
 
-
-
  END SUBROUTINE lsm_check_parameters
 
@@ -1607,24 +1612,25 @@
     LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
 
-    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
-                f1,          & !< resistance correction term 1
-                f2,          & !< resistance correction term 2
-                f3,          & !< resistance correction term 3
-                m_min,       & !< minimum soil moisture
-                e,           & !< water vapour pressure
-                e_s,         & !< water vapour saturation pressure
-                e_s_dt,      & !< derivate of e_s with respect to T
-                tend,        & !< tendency
-                dq_s_dt,     & !< derivate of q_s with respect to T
-                coef_1,      & !< coef. for prognostic equation
-                coef_2,      & !< coef. for prognostic equation
-                f_qsws,      & !< factor for qsws
-                f_qsws_veg,  & !< factor for qsws_veg
-                f_qsws_soil, & !< factor for qsws_soil
-                f_qsws_liq,  & !< factor for qsws_liq
-                f_shf,       & !< factor for shf
-                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
+    REAL(wp) :: c_surface_tmp,  & !< temporary variable for storing the volumetric heat capacity of the surface
+                f1,             & !< resistance correction term 1
+                f2,             & !< resistance correction term 2
+                f3,             & !< resistance correction term 3
+                m_min,          & !< minimum soil moisture
+                e,              & !< water vapour pressure
+                e_s,            & !< water vapour saturation pressure
+                e_s_dt,         & !< derivate of e_s with respect to T
+                tend,           & !< tendency
+                dq_s_dt,        & !< derivate of q_s with respect to T
+                coef_1,         & !< coef. for prognostic equation
+                coef_2,         & !< coef. for prognostic equation
+                f_qsws,         & !< factor for qsws
+                f_qsws_veg,     & !< factor for qsws_veg
+                f_qsws_soil,    & !< factor for qsws_soil
+                f_qsws_liq,     & !< factor for qsws_liq
+                f_shf,          & !< factor for shf
+                lambda_soil,    & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
                 lambda_surface, & !< Current value of lambda_surface (W/m2/K)
-                m_liq_max      !< maxmimum value of the liq. water reservoir
+                m_liq_max,      & !< maxmimum value of the liq. water reservoir
+                ueff              !< limited near-surface wind speed - used for calculation of resistance
 
     TYPE(surf_type_lsm), POINTER ::  surf_t_surface
@@ -1795,14 +1801,18 @@
 !--    to positive values.
        IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
-          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
+          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /                               &
                              ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
        ELSE
-          surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                      &
-                        * ( 11.8_wp + 4.2_wp *                                 &
-                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
-                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
-                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
-                              0.01_wp ) )                                      &
-                           )  - 4.0_wp  )
+!
+!--       Limit wind velocity in order to avoid division by zero.
+!--       The nominator can become <= 0.0 for values z0 < 3*10E-4.
+          ueff        = MAX( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +                    &
+                                    ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +                   &
+                                    ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),                  &
+                             1.0_wp / 4.2_wp * ( 4.0_wp / ( surf%z0(m) * 1000.0_wp ) - 11.8_wp ), &
+                             0.1_wp                                                               &
+                           )
+          surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                                         &
+                                              * ( 11.8_wp + 4.2_wp * ueff )  - 4.0_wp  )
        ENDIF
 !
@@ -2020,5 +2030,5 @@
        ENDIF
 
-       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exner(nzb)
+       surf%pt_surface(m) = surf_t_surface_p%var_1d(m) / exner(nzb)
 
 !
@@ -2034,8 +2044,7 @@
 
        surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
-                                             - surf_t_soil%var_2d(nzb_soil,m) )
+                                      - surf_t_soil%var_2d(nzb_soil,m) )
 
        surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
-
 !
 ! update the 3d field of rad_lw_out array to have consistent output
@@ -2198,5 +2207,4 @@
        CALL calc_q_surface
     ENDIF
-
 !
 !-- Calculate new roughness lengths (for water surfaces only)
@@ -2721,6 +2729,9 @@
                            building_type_f%var(j,i) /= building_type_f%fill )  &
                       THEN
-                         vegetation_type_f%var(j,i)                 = 1 ! bare soil
-                         soil_type_f%var_2d(j,i)                    = 1
+                         vegetation_type_f%var(j,i) = 1 ! bare soil
+                         soil_type_f%var_2d(j,i)    = 1
+
+                         water_type_f%var(j,i)      = water_type_f%fill
+                         pavement_type_f%var(j,i)   = pavement_type_f%fill
 !
 !--                      If surface_fraction is provided in static input,
@@ -2738,7 +2749,7 @@
 !--                Normally proceed with setting surface types.
                    i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
-                                            surf_lsm_v(l)%building_covered(m) )
+                                                   surf_lsm_v(l)%building_covered(m) )
                    j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
-                                            surf_lsm_v(l)%building_covered(m) )
+                                                   surf_lsm_v(l)%building_covered(m) )
                    IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
                       surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
@@ -3821,5 +3832,4 @@
              nzt_pavement = pavement_depth_level
           ENDIF
-
        ENDIF
 !
@@ -4904,8 +4914,6 @@
           DO  m = 1, surf_lsm_h%ns
              IF ( surf_lsm_h%vegetation_surface(m) )  THEN
-                i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,            &
-                                             surf_lsm_v(l)%building_covered(m) )
-                j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff,            &
-                                             surf_lsm_v(l)%building_covered(m) )
+                i = surf_lsm_h%i(m)
+                j = surf_lsm_h%j(m)
                 DO  k = nzb_soil, nzt_soil
                    surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
@@ -5453,6 +5461,5 @@
 !--          true pavement depth is considered
              DO  k = nzb_soil, nzt_soil-1
-                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
-                                        * 0.5_wp
+                surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
              ENDDO
              surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
@@ -5462,8 +5469,9 @@
              tend(:) = 0.0_wp
 
-             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
-                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
-                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) &
-                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
+             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *                          &
+                    ( surf%lambda_h(nzb_soil,m) *                                                  &
+                      ( surf_t_soil%var_2d(nzb_soil+1,m) - surf_t_soil%var_2d(nzb_soil,m) )        &
+                      * ddz_soil_center(nzb_soil)                                                  &
+                    + surf%ghf(m) ) * ddz_soil(nzb_soil)
 
              DO  k = nzb_soil+1, nzt_soil
@@ -5472,5 +5480,5 @@
                      * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
                      * ddz_soil_center(k)                                      &
-                     - surf%lambda_h(k-1,m)                                    &
+                              - surf%lambda_h(k-1,m)                           &
                      * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
                      * ddz_soil_center(k-1)                                    &
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4629)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4630)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! - Bugfix in resistance calculation - avoid potential divisions by zero
+! - Minor formatting adjustment
+!
+! 4602 2020-07-14 14:49:45Z suehring
 ! Add missing initialization of albedo type with values given from static input file
 ! 
@@ -463,4 +467,5 @@
     INTEGER(iwp) ::  wall_category = 2               !< default category for wall surface over pedestrian zone
 
+    REAL(wp)     ::  d_roughness_concrete            !< inverse roughness length of average concrete surface
     REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
 
@@ -3583,4 +3588,7 @@
     ENDIF
 !
+!-- Calculate constant values
+    d_roughness_concrete = 1.0_wp / roughness_concrete
+!
 !-- Flag surface elements belonging to the ground floor level. Therefore, use terrain height array
 !-- from file, if available. This flag is later used to control initialization of surface attributes.
@@ -7539,5 +7547,6 @@
                   q_s,                     &  !< saturation specific humidity
                   rho_lv,                  &  !< frequently used parameter for green layers
-                  tend                      !< tendency
+                  tend,                    &  !< tendency
+                  ueff                        !< limited near-surface wind speed - used for calculation of resistance
 
 
@@ -8046,11 +8055,17 @@
 !--        normal component is obtained by simple linear interpolation. (An alternative would be an
 !--        logarithmic interpolation.) Parameter roughness_concrete (default value = 0.001) is used
-!--        to calculation of roughness relative to concrete
-           surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) / roughness_concrete              &
-                                  * ( 11.8_wp + 4.2_wp                                             &
-                                      * SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +     &
-                                          ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +              &
-                                          ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, 0.01_wp ) )   &
-                                    )  - 4.0_wp  )
+!--        to calculation of roughness relative to concrete. Note, wind velocity is limited
+!--        to avoid division by zero. The nominator can become <= 0.0 for values z0 < 3*10E-4.
+           ueff        = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +                   &
+                                     ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +                   &
+                                     ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),                  &
+                               1.0_wp / 4.2_wp                                                     &
+                               * ( 4.0_wp / ( surf_usm_v(l)%z0(m) * d_roughness_concrete )         &
+                               - 11.8_wp ),                                                        &
+                               0.1_wp                                                              &
+                              )
+
+           surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) * d_roughness_concrete            &
+                                  * ( 11.8_wp + 4.2_wp * ueff )  - 4.0_wp  )
 !
 !--        Limit aerodynamic resistance
