Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4769)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4770)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Consider basal-area density as an additional sink for momentum in the
+! prognostic equations for momentum and SGS TKE
+!
+! 4768 2020-11-02 19:11:23Z suehring
 ! Enable 3D data output also with 64-bit precision
 !
@@ -293,5 +297,6 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !< leaf area density
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
-    
+
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  bad_s                    !< basal-area density on scalar-grid
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf               !< cumulative lai for heatflux calc.
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
@@ -546,5 +551,5 @@
           unit = 'K s-1'
 
-       CASE ( 'pcm_lad' )
+       CASE ( 'pcm_bad', 'pcm_lad' )
           unit = 'm2 m-3'
 
@@ -870,4 +875,15 @@
           ENDIF
 
+       CASE ( 'pcm_bad' )
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = MAX( 1, nzb_do ), MIN( pch_index_ji(j,i), nzt_do )
+                      local_pf(i,j,k) = bad_s(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
        CASE DEFAULT
           found = .FALSE.
@@ -898,5 +914,5 @@
      SELECT CASE ( TRIM( var ) )
 
-        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
+        CASE ( 'pcm_heatrate', 'pcm_bad', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate')
            grid_x = 'x'
            grid_y = 'y'
@@ -1039,4 +1055,5 @@
 
        LOGICAL      ::  lad_on_top = .FALSE.  !< dummy flag to indicate that LAD is defined on a building roof
+       LOGICAL      ::  bad_on_top = .FALSE.  !< dummy flag to indicate that BAD is defined on a building roof
 
        REAL(wp) ::  canopy_height   !< canopy height for lad-profile construction
@@ -1142,7 +1159,9 @@
 
 !
-!--    Allocate 3D-array for the leaf area density (lad_s).
+!--    Allocate 3D-array for the leaf-area density (lad_s) as well as 
+!--    for basal-area densitiy (bad_s). Note, by default bad_s is zero.
        ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
+       ALLOCATE( bad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       bad_s = 0.0_wp
 !
 !--    Initialization of the canopy coverage in the model domain:
@@ -1303,10 +1322,8 @@
 #endif
                 IF ( lad_on_top )  THEN
-                   WRITE( message_string, * )                                  &
-                                        'Resolved plant-canopy is ' //         &
-                                        'defined on top of an artificially '// &
-                                        'created building grid point(s) ' //   &
-                                        '(emerged from the filtering) - ' //   &
-                                        'LAD profile is omitted at this / ' // &
+                   WRITE( message_string, * )                                                      &
+                                        'Resolved plant-canopy is defined on top of an ' //        &
+                                        'artificially created building grid point(s) '//           &
+                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
                                         'these grid point(s).'
                    CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
@@ -1321,4 +1338,47 @@
              ELSE
                 CALL pcm_read_plant_canopy_3d
+             ENDIF
+!
+!--          Initialize LAD with data from file. If LAD is given in NetCDF file,
+!--          use these values, else take LAD profiles from ASCII file. 
+!--          Please note, in NetCDF file LAD is only given up to the maximum 
+!--          canopy top, indicated by basal_area_density_f%nz.  
+             bad_s = 0.0_wp
+             IF ( basal_area_density_f%from_file )  THEN
+!
+!--             Set also pch_index, used to be the upper bound of the vertical 
+!--             loops. Therefore, use the global top of the canopy layer. 
+                pch_index = basal_area_density_f%nz - 1
+
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = 0, basal_area_density_f%nz - 1
+                         IF ( basal_area_density_f%var(k,j,i) /= basal_area_density_f%fill )       &
+                            bad_s(k,j,i) = basal_area_density_f%var(k,j,i)
+                      ENDDO
+!
+!--                   Check if resolved vegetation is mapped onto buildings.
+!--                   Please see comment for leaf_area density
+                      IF ( ANY( bad_s(:,j,i) /= 0.0_wp )  .AND.                &
+                           ANY( BTEST( wall_flags_total_0(:,j,i), 6 ) ) .AND.  &
+                           ANY( BTEST( wall_flags_total_0(:,j,i), 4 ) ) )  THEN
+                         bad_s(:,j,i) = 0.0_wp
+                         bad_on_top   = .TRUE.
+                      ENDIF
+                   ENDDO
+                ENDDO
+#if defined( __parallel )
+               CALL MPI_ALLREDUCE( MPI_IN_PLACE, bad_on_top, 1, MPI_LOGICAL,  &
+                                   MPI_LOR, comm2d, ierr)
+#endif
+                IF ( bad_on_top )  THEN
+                   WRITE( message_string, * )                                                      &
+                                        'Resolved plant-canopy is defined on top of an ' //        &
+                                        'artificially created building grid point(s) '//           &
+                                        'the filtering) - LAD/BAD profile is omitted at this / ' //&
+                                        'these grid point(s).'
+                   CALL message( 'pcm_init', 'PA0313', 0, 0, 0, 6, 0 )
+                ENDIF
+                CALL exchange_horiz( bad_s, nbgp )
              ENDIF
 
@@ -1355,5 +1415,5 @@
           DO  j = nysg, nyng
              DO  k = 0, pch_index
-                IF ( lad_s(k,j,i) /= 0 )  pch_index_ji(j,i) = k
+                IF ( lad_s(k,j,i) /= 0.0_wp  .OR.  bad_s(k,j,i) /= 0.0_wp )  pch_index_ji(j,i) = k
              ENDDO
 !
@@ -1862,4 +1922,5 @@
        LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
 
+       REAL(wp) ::  bad_local !< local bad value
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
        REAL(wp) ::  lad_local !< local lad value
@@ -1906,8 +1967,13 @@
 !--                   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. 
+!--                   on top of a roof where (usually) is no LAD is defined. The same is
+!--                   also valid for bad_s.
                       lad_local = lad_s(kk,j,i)
-                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  &
+                      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)
+
+                      bad_local = bad_s(kk,j,i)
+                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp                    &
+                           .AND.  .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
 !
 !--                   In order to avoid that LAD is mistakenly considered at right-
@@ -1915,7 +1981,9 @@
 !--                   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  &
+!--                   point (j,i-1). The same is also valid for bad_s.
+                      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)
+                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp                    &
+                           .AND.  building_edge_e )  bad_local = bad_s(kk,j,i-1)
 
                       pre_tend = 0.0_wp
@@ -1924,5 +1992,5 @@
 !--                   Calculate preliminary value (pre_tend) of the tendency
                       pre_tend = - canopy_drag_coeff *                         &
-                                   lad_local *                                 &
+                                   ( lad_local + bad_local ) *                 &
                                    SQRT( u(k,j,i)**2 +                         &
                                          ( 0.25_wp * ( v(k,j,i-1) +            &
@@ -1987,7 +2055,12 @@
 !--                   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. 
+!--                   The same is also valid for bad_s.
                       lad_local = lad_s(kk,j,i)
-                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
+                      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)
+
+                      bad_local = bad_s(kk,j,i)
+                      IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp                    &
+                       .AND.  .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
 !
 !--                   In order to avoid that LAD is mistakenly considered at right-
@@ -1995,7 +2068,10 @@
 !--                   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  &
+!--                   point (j,i-1). The same is also valid for bad_s.
+                      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)
+
+                      IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp                    &
+                       .AND.  building_edge_n )  bad_local = bad_s(kk,j-1,i)
 
                       pre_tend = 0.0_wp
@@ -2004,5 +2080,5 @@
 !--                   Calculate preliminary value (pre_tend) of the tendency
                       pre_tend = - canopy_drag_coeff *                         &
-                                   lad_local *                                 &
+                                   ( lad_local + bad_local ) *                 &
                                    SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
                                                        u(k,j-1,i+1) +          &
@@ -2054,6 +2130,8 @@
 !--                   Calculate preliminary value (pre_tend) of the tendency
                       pre_tend = - canopy_drag_coeff *                         &
-                                   (0.5_wp *                                   &
-                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
+                                   ( 0.5_wp *                                  &
+                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +    &
+                                     0.5_wp *                                  &
+                                      ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *  &
                                    SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
                                                        u(k,j,i+1) +            &
@@ -2161,5 +2239,5 @@
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        2.0_wp * canopy_drag_coeff *            &
-                                       lad_s(kk,j,i) *                         &
+                                       ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *     &
                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
                                                           u(k,j,i+1) )         &
@@ -2253,4 +2331,5 @@
        LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
 
+       REAL(wp) ::  bad_local !< local lad value
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
        REAL(wp) ::  lad_local !< local lad value
@@ -2296,7 +2375,12 @@
 !--             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. 
+!--             The same is also valid for bad_s.
                 lad_local = lad_s(kk,j,i)
-                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  .AND. &
+                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)
+
+                bad_local = bad_s(kk,j,i)
+                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j,i-1) > 0.0_wp  .AND.                   &
+                     .NOT. building_edge_w )  bad_local = bad_s(kk,j,i-1)
 !
 !--             In order to avoid that LAD is mistakenly considered at right-
@@ -2304,7 +2388,9 @@
 !--             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. &
+!--             point (j,i-1). The same is also valid for bad_s.
+                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)
+                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j,i-1) == 0.0_wp  .AND.                   &
+                     building_edge_e )  bad_local = bad_s(kk,j,i-1)
 
                 pre_tend = 0.0_wp
@@ -2313,5 +2399,5 @@
 !--             Calculate preliminary value (pre_tend) of the tendency
                 pre_tend = - canopy_drag_coeff *                               &
-                             lad_local *                                       &   
+                             ( lad_local + bad_local ) *                       &
                              SQRT( u(k,j,i)**2 +                               &
                                    ( 0.25_wp * ( v(k,j,i-1)  +                 &
@@ -2374,7 +2460,12 @@
 !--             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. 
+!--             The same is also valid for bad_s.
                 lad_local = lad_s(kk,j,i)
-                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp  .AND. &
+                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)
+
+                bad_local = bad_s(kk,j,i)
+                IF ( bad_local == 0.0_wp  .AND.  bad_s(kk,j-1,i) > 0.0_wp  .AND.                   &
+                     .NOT. building_edge_s )  bad_local = bad_s(kk,j-1,i)
 !
 !--             In order to avoid that LAD is mistakenly considered at right-
@@ -2382,7 +2473,9 @@
 !--             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. &
+!--             point (j,i-1). The same is also valid for bad_s.
+                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)
+                IF ( bad_local > 0.0_wp  .AND.  bad_s(kk,j-1,i) == 0.0_wp  .AND.                   &
+                     building_edge_n )  bad_local = bad_s(kk,j-1,i)
 
                 pre_tend = 0.0_wp
@@ -2391,5 +2484,5 @@
 !--             Calculate preliminary value (pre_tend) of the tendency
                 pre_tend = - canopy_drag_coeff *                               &
-                             lad_local *                                       &
+                             ( lad_local + bad_local ) *                       &
                              SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
                                                  u(k,j-1,i+1) +                &
@@ -2437,7 +2530,9 @@
 !--             Calculate preliminary value (pre_tend) of the tendency
                 pre_tend = - canopy_drag_coeff *                               &
-                             (0.5_wp *                                         &
-                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
-                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &  
+                             ( 0.5_wp *                                        &
+                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) ) +          &
+                               0.5_wp *                                        &
+                                ( bad_s(kk+1,j,i) + bad_s(kk,j,i) ) ) *        &
+                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &
                                                  u(k,j,i+1)  +                 &
                                                  u(k+1,j,i)  +                 &
@@ -2527,5 +2622,5 @@
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  2.0_wp * canopy_drag_coeff *                  &
-                                 lad_s(kk,j,i) *                               &
+                                 ( lad_s(kk,j,i) + bad_s(kk,j,i) ) *           &
                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
                                                     u(k,j,i+1) )               &
@@ -2540,5 +2635,4 @@
                                  e(k,j,i)
              ENDDO
-             
 !
 !--       scalar concentration 
