Index: /palm/trunk/SOURCE/init_grid.f90
===================================================================
--- /palm/trunk/SOURCE/init_grid.f90	(revision 4313)
+++ /palm/trunk/SOURCE/init_grid.f90	(revision 4314)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Set additional topography flag 4 to mark topography grid points emerged
+! from the filtering process.
+! 
+! 4294 2019-11-13 18:34:16Z suehring
 ! Bugfix, always set bit 5 and 6 of wall_flags, indicating terrain- and 
 ! building surfaces in all  cases, in order to enable terrain-following output 
@@ -1600,8 +1604,8 @@
                       num_hole_l     = num_hole_l + 1
 !
-!--                   Clear flag 0 and set special flag ( bit 3) to indicate 
+!--                   Clear flag 0 and set special flag ( bit 4) to indicate 
 !--                   that new topography point is a result of filtering process.
                       topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
-                      topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 3 )
+                      topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 4 )
 !
 !--                   If filled grid point is occupied by a building, classify
@@ -2180,5 +2184,5 @@
           CALL exchange_horiz_int( topo, nys, nyn, nxl, nxr, nzt, nbgp )
 !
-!--       Set lateral boundary conditions for topography on all ghost layers          
+!--       Set lateral boundary conditions for topography on all ghost layers
           IF ( .NOT. bc_ns_cyc )  THEN
              IF ( nys == 0  )  THEN
@@ -2601,4 +2605,14 @@
     ENDDO
 !
+!-- Set flag 4, indicating new topography grid points due to filtering.
+    DO i = nxl, nxr
+       DO j = nys, nyn
+          DO k = nzb, nzt+1
+             IF ( BTEST( topo(k,j,i), 4 ) )                                    &
+                wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
+          ENDDO
+       ENDDO
+    ENDDO
+!
 !-- Exchange ghost points for wall flags
     CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4313)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4314)
@@ -27,4 +27,11 @@
 ! -----------------
 ! $Id$
+! - Bugfix, plant canopy was still considered at building edges on for the u- 
+!   and v-component.
+! - Relax restriction of LAD on building tops. LAD is only omitted at 
+!   locations where building grid points emerged artificially by the 
+!   topography filtering.
+! 
+! 4309 2019-11-26 18:49:59Z suehring
 ! Typo
 ! 
@@ -114,5 +121,5 @@
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
-               nz, nzb, nzt, topo_top_ind
+               nz, nzb, nzt, topo_top_ind, wall_flags_0
 
     USE kinds
@@ -938,7 +945,4 @@
            ONLY: message_string, ocean_mode
 
-       USE indices,                                                            &
-           ONLY:  wall_flags_0
-
        USE netcdf_data_input_mod,                                              &
            ONLY:  leaf_area_density_f
@@ -1108,17 +1112,27 @@
 !
 !--                   Check if resolved vegetation is mapped onto buildings.
-!--                   Even if the user defines LAD only for non-building covered
-!--                   grid points, the filter process of topography may cause
-!--                   that LAD overlaps with buildings. This case, an 
-!--                   informative message is given and the LAD field is set to
-!--                   zero at these 
+!--                   In general, this is allowed and also meaningful, e.g. 
+!--                   when trees carry across roofs. However, due to the 
+!--                   topography filtering, new building grid points can emerge
+!--                   at location where also plant canopy is defined. As a 
+!--                   result, plant canopy is mapped on top of roofs, with
+!--                   siginficant impact on the downstream flow field and the
+!--                   nearby surface radiation. In order to avoid that 
+!--                   plant canopy is mistakenly mapped onto building roofs, 
+!--                   check for building grid points (bit 6) that emerge from 
+!--                   the filtering (bit 4) and set LAD to zero at these
+!--                   artificially  created building grid points. This case, 
+!--                   an informative message is given.
                       IF ( ANY( lad_s(:,j,i) /= 0.0_wp )  .AND.                &
-                           ANY( BTEST( wall_flags_0(:,j,i), 6 ) ) )  THEN
+                           ANY( BTEST( wall_flags_0(:,j,i), 6 ) ) .AND.        &
+                           ANY( BTEST( wall_flags_0(:,j,i), 4 ) ) )  THEN
                          lad_s(:,j,i) = 0.0_wp
                          WRITE( message_string, * )                            &
-                                          'Resolved plant-canopy is ' //       &
-                                          'defined on top of a building ' //   &
-                                          '- lad is omitted at this grid '   //&
-                                          'point: (i,j) = ', i, j
+                                        'Resolved plant-canopy is ' //         &
+                                        'defined on top of an artificially '// &
+                                        'created building grid point ' //      &
+                                        '(emerged from the filering) - ' //    &
+                                        'LAD profile is omitted at this ' //   &
+                                        'grid point: (i,j) = ', i, j
                          CALL message( 'pcm_init', 'PA0313', 0, 0, myid, 6, 0 )
                       ENDIF
@@ -1918,4 +1932,9 @@
        INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
+       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
+       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
+       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
+       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
+
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
        REAL(wp) ::  lad_local !< local lad value
@@ -1935,4 +1954,13 @@
           CASE ( 1 )
 !
+!--          Set control flags indicating east- and westward-orientated 
+!--          building edges. Note, building_egde_w is set from the perspective
+!--          of the potential rooftop grid point, while building_edge_e is
+!--          is set from the perspective of the non-building grid point.
+             building_edge_w = ANY( BTEST( wall_flags_0(:,j,i),   6 ) )  .AND. &
+                         .NOT. ANY( BTEST( wall_flags_0(:,j,i-1), 6 ) )
+             building_edge_e = ANY( BTEST( wall_flags_0(:,j,i-1), 6 ) )  .AND. &
+                         .NOT. ANY( BTEST( wall_flags_0(:,j,i),   6 ) )
+!
 !--          Determine topography-top index on u-grid
              k_wall = topo_top_ind(j,i,1)
@@ -1948,9 +1976,20 @@
 !--             than inside of the canopy.
 !--             For the same reason, the lad at the rightmost(i+1)canopy 
-!--             boundary on the u-grid equals lad_s(k,j,i).
+!--             boundary on the u-grid equals lad_s(k,j,i), which is considered
+!--             in the next if-statement. Note, at left-sided building edges 
+!--             this is not applied, here the LAD is equals the LAD at grid
+!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped 
+!--             on top of a roof where (usually) is no LAD is defined. 
                 lad_local = lad_s(kk,j,i)
-                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
-                   lad_local = lad_s(kk,j,i-1)
-                ENDIF
+                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  .AND. &
+                     .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
+!
+!--             In order to avoid that LAD is mistakenly considered at right-
+!--             sided building edges (here the topography-top index for the
+!--             u-component at index j,i is still on the building while the
+!--             topography top for the scalar isn't), LAD is taken from grid
+!--             point (j,i-1). 
+                IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  .AND. &
+                     building_edge_e )  lad_local = lad_s(kk,j,i-1)
 
                 pre_tend = 0.0_wp
@@ -1995,7 +2034,15 @@
           CASE ( 2 )
 !
+!--          Set control flags indicating north- and southward-orientated 
+!--          building edges. Note, building_egde_s is set from the perspective
+!--          of the potential rooftop grid point, while building_edge_n is
+!--          is set from the perspective of the non-building grid point.
+             building_edge_s = ANY( BTEST( wall_flags_0(:,j,i),   6 ) )  .AND. &
+                         .NOT. ANY( BTEST( wall_flags_0(:,j-1,i), 6 ) )
+             building_edge_n = ANY( BTEST( wall_flags_0(:,j-1,i), 6 ) )  .AND. &
+                         .NOT. ANY( BTEST( wall_flags_0(:,j,i),   6 ) )
+!
 !--          Determine topography-top index on v-grid
              k_wall = topo_top_ind(j,i,2)
-
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
 
@@ -2008,9 +2055,20 @@
 !--             than inside of the canopy.
 !--             For the same reason, the lad at the northmost(j+1)canopy 
-!--             boundary on the v-grid equals lad_s(k,j,i).
+!--             boundary on the v-grid equals lad_s(k,j,i), which is considered
+!--             in the next if-statement. Note, at left-sided building edges 
+!--             this is not applied, here the LAD is equals the LAD at grid
+!--             point (k,j,i), in order to avoid that LAD is mistakenly mapped 
+!--             on top of a roof where (usually) is no LAD is defined. 
                 lad_local = lad_s(kk,j,i)
-                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
-                   lad_local = lad_s(kk,j-1,i)
-                ENDIF
+                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp  .AND. &
+                     .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
+!
+!--             In order to avoid that LAD is mistakenly considered at right-
+!--             sided building edges (here the topography-top index for the
+!--             u-component at index j,i is still on the building while the
+!--             topography top for the scalar isn't), LAD is taken from grid
+!--             point (j,i-1). 
+                IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  .AND. &
+                     building_edge_n )  lad_local = lad_s(kk,j-1,i)
 
                 pre_tend = 0.0_wp
