Index: /palm/trunk/SOURCE/pmc_handle_communicator_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_handle_communicator_mod.f90	(revision 4770)
+++ /palm/trunk/SOURCE/pmc_handle_communicator_mod.f90	(revision 4771)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Canopy-restricted anterpolation introduced. New namelist parameter anterpolation_starting_height 
+! introduced for controlling canopy-restricted anterpolation.
+!
+! 4650 2020-08-25 14:35:50Z raasch
 ! bugfix for r4649
 !
@@ -138,5 +142,5 @@
 !--------------------------------------------------------------------------------------------------!
  SUBROUTINE pmc_init_model( comm, nesting_datatransfer_mode, nesting_mode,                         &
-                            anterpolation_buffer_width, pmc_status )
+                            anterpolation_buffer_width, anterpolation_starting_height, pmc_status )
 
     USE control_parameters,                                                                        &
@@ -155,4 +159,6 @@
     INTEGER, INTENT(INOUT) ::  pmc_status                  !<
 
+    REAL(wp), INTENT(INOUT) ::  anterpolation_starting_height  !< steering parameter for canopy restricted anterpolation 
+    
     INTEGER ::  childcount     !<
     INTEGER ::  i              !<
@@ -180,5 +186,6 @@
 
        CALL read_coupling_layout( nesting_datatransfer_mode, nesting_mode,                         &
-                                  anterpolation_buffer_width, pmc_status )
+            anterpolation_buffer_width, anterpolation_starting_height,                             &
+            pmc_status )
 
        IF ( pmc_status /= pmc_no_namelist_found  .AND.                                             &
@@ -252,4 +259,5 @@
                     MPI_COMM_WORLD, istat )
     CALL MPI_BCAST( anterpolation_buffer_width, 1, MPI_INT, 0, MPI_COMM_WORLD, istat )
+    CALL MPI_BCAST( anterpolation_starting_height, 1, MPI_REAL, 0, MPI_COMM_WORLD, istat )
 !
 !-- Assign global MPI processes to individual models by setting the couple id
@@ -422,5 +430,6 @@
 !--------------------------------------------------------------------------------------------------!
  SUBROUTINE read_coupling_layout( nesting_datatransfer_mode, nesting_mode,                         &
-                                  anterpolation_buffer_width, pmc_status )
+                                  anterpolation_buffer_width, anterpolation_starting_height,       &
+                                  pmc_status )
 
     IMPLICIT NONE
@@ -431,4 +440,6 @@
     INTEGER, INTENT(INOUT)      ::  anterpolation_buffer_width  !< Boundary buffer width for anterpolation
     INTEGER(iwp), INTENT(INOUT) ::  pmc_status                  !<
+
+    REAL(wp), INTENT(INOUT) ::  anterpolation_starting_height   !< steering parameter for canopy restricted anterpolation
 
     INTEGER(iwp) ::  bad_llcorner  !<
@@ -441,5 +452,6 @@
                                    nesting_datatransfer_mode,                                      &
                                    nesting_mode,                                                   &
-                                   anterpolation_buffer_width
+                                   anterpolation_buffer_width,                                     &
+                                   anterpolation_starting_height
 
 !
Index: /palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4770)
+++ /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 4771)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Canopy-restricted anterpolation introduced. New namelist parameter anterpolation_starting_height
+! introduced for controlling canopy-restricted anterpolation.
+!
+! 4745 2020-10-15 16:37:13Z suehring
 ! Adjustement of face-area calculation to 3D-topographies
 !
@@ -453,4 +457,6 @@
     INTEGER(iwp), SAVE ::  anterpolation_buffer_width = 2  !< Boundary buffer width for anterpolation
 
+    REAL(wp), SAVE ::  anterpolation_starting_height = 9999999.9_wp  !< steering parameter for canopy restricted anterpolation
+
     LOGICAL, SAVE ::  nested_run = .FALSE.        !< general switch
     LOGICAL, SAVE ::  rans_mode_parent = .FALSE.  !< mode of parent model (.F. - LES mode, .T. - RANS mode)
@@ -463,4 +469,5 @@
 !
 !-- Children's parent-grid arrays
+    INTEGER(iwp), SAVE, DIMENSION(:,:), ALLOCATABLE :: kpb_anterp  !< Lower limit of kp in anterpolation 
     INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC ::  parent_bound  !< subdomain index bounds for children's parent-grid arrays
 
@@ -702,5 +709,5 @@
 
     CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,                      &
-                         anterpolation_buffer_width, pmc_status )
+                         anterpolation_buffer_width, anterpolation_starting_height, pmc_status )
 
     IF ( pmc_status == pmc_no_namelist_found )  THEN
@@ -1507,5 +1514,7 @@
 !--    Compute surface areas of the nest-boundary faces
        CALL pmci_compute_face_areas
-
+!       
+!--    Compute anterpolation lower kp-index bounds if necessary
+       CALL pmci_compute_kpb_anterp 
     ENDIF
 
@@ -2438,4 +2447,155 @@
 
  END SUBROUTINE pmci_compute_face_areas
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Define the anterpolation starting kp-indices for all jp,ip based on the obstacle-canopy 
+!> topograhy for canopy-restricted anterpolation. Note that this is based on the child terrain
+!> topography information since it is difficult to access the parent topography information 
+!> from the child. This means that these topographies are assumed to be close to each other. 
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE pmci_compute_kpb_anterp
+
+    INTEGER(iwp) ::  ic                     !< Child-grid index in the x-direction
+    INTEGER(iwp) ::  ip                     !< Parent-grid index in the x-direction
+    INTEGER(iwp) ::  jc                     !< Child-grid index in the y-direction
+    INTEGER(iwp) ::  jp                     !< Parent-grid index in the y-direction
+    INTEGER(iwp) ::  terrain_surf_k_child   !< Local terrain height k index in the child grid
+    INTEGER(iwp) ::  terrain_surf_k_parent  !< Local terrain height k index in the parent grid
+
+    
+    ALLOCATE( kpb_anterp(jps:jpn,ipl:ipr) )
+    kpb_anterp(jps:jpn,ipl:ipr) = 0
+   
+    IF ( nesting_mode /= 'one-way'  .AND.  topography /= 'flat' )  THEN
+!      
+!--    Check if the value was given by user or not
+       IF ( anterpolation_starting_height > 9000000.0_wp )  THEN
+!         
+!--       No value was given by user, compute a default value first.
+          CALL default_anterpolation_starting_height( 99 )
+       ENDIF
+
+       DO  ip = ipl, ipr
+          DO  jp = jps, jpn            
+             terrain_surf_k_parent = 0
+             DO  ic = iflo(ip), ifuo(ip)
+                DO  jc = jflo(jp), jfuo(jp)
+                   terrain_surf_k_child = MINLOC(                                                   &
+                        MERGE( 1, 0, BTEST( wall_flags_total_0(:,jc,ic), 5 ) ), DIM = 1 ) - 1
+                   terrain_surf_k_parent = MAX( terrain_surf_k_child, terrain_surf_k_parent )
+                ENDDO
+             ENDDO
+             terrain_surf_k_parent = NINT( terrain_surf_k_parent / REAL( kgsr, wp ) )
+             kpb_anterp(jp,ip) = terrain_surf_k_parent + NINT( anterpolation_starting_height / pg%dz )
+          ENDDO
+       ENDDO
+
+    ENDIF
+
+ END SUBROUTINE pmci_compute_kpb_anterp
+
+
+!--------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Compute a default value for the anterpolation_starting_height for the canopy-restricted 
+!> anterpolation. The default value is based on a given percentile of the child building height 
+!> distribution. This operation is based on the child-grid information.
+!--------------------------------------------------------------------------------------------------!
+ SUBROUTINE default_anterpolation_starting_height( percentile_level ) 
+   
+    INTEGER(iwp), INTENT(IN) ::  percentile_level  !< Selected percentile level (1,...,99)
+
+    INTEGER(iwp) ::  ic                            !< Child-grid index in the x-direction
+    INTEGER(iwp) ::  ierr                          !< MPI-error code
+    INTEGER(iwp) ::  im                            !< Grid index in the x-direction of the currently
+                                                   !< globally highest building node
+    INTEGER(iwp) ::  jc                            !< Child-grid index in the y-direction
+    INTEGER(iwp), DIMENSION(2) :: ji_maxloc        !< jc and ic indices of the local maximum height
+    INTEGER(iwp) ::  jm                            !< Grid index in the y-direction of the currently
+                                                   !< globally highest building node
+    INTEGER(iwp) ::  num_building_nodes            !< Global number of jc,ic-nodes under buildings
+    INTEGER(iwp) ::  num_building_nodes_l          !< Local number of jc,ic-nodes under buildings
+    INTEGER(iwp) ::  num_highest_nodes             !< The number of highest building nodes above the
+                                                   !< percentile value
+    INTEGER(iwp) ::  n                             !< Running index for the sorting loop 
+    INTEGER(iwp) ::  pid_highest                   !< Process id of the process where the currently
+                                                   !< highest building node was found.
+    INTEGER(iwp) ::  terrain_surf_k                !< Local terrain height k index in the child grid
+                                                  
+    REAL(wp) ::  global_currently_highest          !< Globally highest building node currently in the
+                                                   !< search loop
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  building_height  !< Temporary array for buiding heights
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  max_height         !< Maximum buiding heights of each 
+                                                               !< process 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  max_height_l       !< max_height_l( myid ) is the local 
+                                                               !< maximum buiding height of the 
+                                                               !< current process while the other 
+                                                               !< elements are zeroes
+    
+    ALLOCATE( building_height(nys:nyn,nxl:nxr) )
+!
+!-- Find the building heights and temporarily store them in building_height(nys:nyn,nxl:nxr).    
+!-- Also find the number of nodes occupied by buildings
+    num_building_nodes_l = 0
+    num_building_nodes   = 0
+    DO  ic = nxl, nxr
+       DO  jc = nys, nyn
+          terrain_surf_k  = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jc,ic), 5 ) ),         &
+                                    DIM = 1 ) - 1
+          IF ( topo_top_ind(jc,ic,5) > terrain_surf_k )  THEN                      
+             building_height(jc,ic) = dz(1) * ( topo_top_ind(jc,ic,5) - terrain_surf_k + 1 )
+             num_building_nodes_l = num_building_nodes_l + 1
+          ENDIF
+       ENDDO
+    ENDDO
+    
+    CALL MPI_ALLREDUCE( num_building_nodes_l, num_building_nodes, 1, MPI_INT, MPI_SUM, comm2d, ierr )
+    
+    num_highest_nodes = MAX( 1, NINT( num_building_nodes * ( 100 - percentile_level ) / 100.0_wp ) )
+!    
+!-- Find global percentile value
+    ALLOCATE( max_height_l(0:numprocs-1) )
+    ALLOCATE( max_height(0:numprocs-1) )
+    max_height_l = 0.0_wp
+!
+!-- Search loop for the num_highest_nodes highest nodes globally.
+    DO  n = 1, num_highest_nodes
+       max_height = 0.0_wp
+       max_height_l(myid) = MAXVAL( building_height )
+       CALL MPI_ALLREDUCE( max_height_l, max_height, numprocs, MPI_REAL, MPI_SUM, comm2d, ierr )
+       global_currently_highest = MAXVAL( max_height )
+       pid_highest = MAXLOC( max_height, DIM = 1 ) - 1
+!       
+!--    Set the globally currently highest height to zero in order to find the next 
+!--    highest node in the next loop cycle.
+       IF ( pid_highest == myid )  THEN
+          ji_maxloc = MAXLOC( building_height )
+          jm = ji_maxloc(1) - 1 + nys 
+          im = ji_maxloc(2) - 1 + nxl
+          building_height(jm,im) = 0.0_wp
+       ENDIF
+    ENDDO
+!
+!-- Now after the search-loop, global_currently_highest is the desired percentile 
+!-- value. Add two parent-grid dz in order to create appropriate clearance above the
+!-- roofs in cases of evenly distributed building height such as simple 
+!-- cuboid-array test cases.
+    anterpolation_starting_height = global_currently_highest + 2.0_wp * pg%dz
+
+    WRITE( 9, * )
+    WRITE( 9, "('A default value is defined and set for anterpolation_starting_height = ',          &
+                f4.1, ' m.' )" )  anterpolation_starting_height
+    WRITE( 9, * )
+
+    DEALLOCATE( max_height_l )
+    DEALLOCATE( max_height )
+    DEALLOCATE( building_height )
+    
+ END SUBROUTINE default_anterpolation_starting_height
+
 #endif
 
@@ -4815,5 +4975,4 @@
     INTEGER(iwp) ::  kc              !< Running index z-direction - child grid
     INTEGER(iwp) ::  kp              !< Running index z-direction - parent grid
-    INTEGER(iwp) ::  kpb_anterp = 0  !< Bottom boundary index for anterpolation along z
     INTEGER(iwp) ::  kpt_anterp      !< Top boundary index for anterpolation along z
     INTEGER(iwp) ::  var_flag        !< bit number used to flag topography on respective grid
@@ -4834,5 +4993,7 @@
     jps_anterp = jps
     jpn_anterp = jpn
-    kpb_anterp = 0
+!
+!-- kpb_anterp is a function of jp and ip and it is set in the initialization phase in 
+!-- pmci_compute_kpb_anterp.
     kpt_anterp = kct - 1 - anterpolation_buffer_width
 
@@ -4862,6 +5023,8 @@
        DO  jp = jps_anterp, jpn_anterp
 !
-!--       For simplicity anterpolate within buildings and under elevated terrain too
-          DO  kp = kpb_anterp, kpt_anterp
+!--       If the user has set anterpolation_starting_height less than the canopy height, the 
+!--       anterpolation is made also within buildings for simplicity, and even under elevated 
+!--       terrain if anterpolation_starting_height is set smaller than terrain height.
+          DO  kp = kpb_anterp(jp,ip), kpt_anterp
              cellsum = 0.0_wp
              DO  ic = ifl(ip), ifu(ip)
