Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2601)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 2602)
@@ -26,5 +26,7 @@
 ! -----------------
 ! $Id$
-! Some cleanup and commenting improvements only.
+! Index-limit related bug (occurred with nesting_mode='vertical') fixed in
+! pmci_interp_tril_t. Check for too high nest domain added in pmci_setup_parent.   
+! Some cleaning up made. 
 ! 
 ! 2582 2017-10-26 13:19:46Z hellstea
@@ -159,5 +161,5 @@
 ! only the total number of PEs is given for the domains, npe_x and npe_y
 ! replaced by npe_total, two unused elements removed from array
-! define_coarse_grid_real,
+! parent_grid_info_real,
 ! array management changed from linked list to sequential loop
 !
@@ -258,10 +260,8 @@
 
     PRIVATE
-
 !
 !-- Constants
     INTEGER(iwp), PARAMETER ::  child_to_parent = 2   !:
     INTEGER(iwp), PARAMETER ::  parent_to_child = 1   !:
-
 !
 !-- Coupler setup
@@ -271,5 +271,4 @@
     INTEGER(iwp), SAVE      ::  cpl_npe_total          !:
     INTEGER(iwp), SAVE      ::  cpl_parent_id          !:
-
 !
 !-- Control parameters, will be made input parameters later
@@ -287,5 +286,4 @@
     REAL(wp), SAVE ::  anterp_relax_length_n = -1.0_wp   !:
     REAL(wp), SAVE ::  anterp_relax_length_t = -1.0_wp   !:
-
 !
 !-- Geometry
@@ -295,5 +293,4 @@
     REAL(wp), SAVE                            ::  lower_left_coord_x   !:
     REAL(wp), SAVE                            ::  lower_left_coord_y   !:
-
 !
 !-- Child coarse data arrays
@@ -321,5 +318,4 @@
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !:
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc   !:
-
 !
 !-- Child interpolation coefficients and child-array indices to be 
@@ -343,5 +339,4 @@
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
-
 !
 !-- Child index arrays and log-ratio arrays for the log-law near-wall
@@ -372,10 +367,8 @@
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
-
 !
 !-- Upper bounds for k in anterpolation.
     INTEGER(iwp), SAVE ::  kctu   !:
     INTEGER(iwp), SAVE ::  kctw   !:
-
 !
 !-- Upper bound for k in log-law correction in interpolation.
@@ -384,5 +377,4 @@
     INTEGER(iwp), SAVE ::  nzt_topo_nestbc_r   !:
     INTEGER(iwp), SAVE ::  nzt_topo_nestbc_s   !:
-
 !
 !-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
@@ -391,5 +383,4 @@
     INTEGER(iwp), SAVE ::  nhls   !:
     INTEGER(iwp), SAVE ::  nhln   !:
-
 !
 !-- Spatial under-relaxation coefficients for anterpolation.
@@ -397,5 +388,4 @@
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
-
 !
 !-- Child-array indices to be precomputed and stored for anterpolation.
@@ -412,5 +402,4 @@
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
-
 !
 !-- Number of fine-grid nodes inside coarse-grid ij-faces
@@ -422,6 +411,6 @@
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_s         !:
     
-    INTEGER(iwp), DIMENSION(3)          ::  define_coarse_grid_int    !:
-    REAL(wp), DIMENSION(7)              ::  define_coarse_grid_real   !:
+    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !:
+    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !:
 
     TYPE coarsegrid_def
@@ -535,5 +524,4 @@
 
     ENDIF
-
 !
 !-- Check steering parameter values
@@ -554,9 +542,7 @@
        CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
     ENDIF
-
 !
 !-- Set the general steering switch which tells PALM that its a nested run
     nested_run = .TRUE.
-
 !
 !-- Get some variables required by the pmc-interface (and in some cases in the
@@ -650,4 +636,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_ys   !:
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_yn   !:
+    REAL(wp) ::  cl_height        !: 
     REAL(wp) ::  dx_cl            !:
     REAL(wp) ::  dy_cl            !:
@@ -663,10 +650,8 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !:
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !:
-    
 
 !
 !   Initialize the pmc parent
     CALL pmc_parentinit
-
 !
 !-- Corners of all children of the present parent
@@ -677,5 +662,4 @@
        ALLOCATE( ch_yn(1:SIZE( pmc_parent_for_child ) - 1) )
     ENDIF
-
 !
 !-- Get coordinates from all children
@@ -692,7 +676,6 @@
           dx_cl = val(4)
           dy_cl = val(5)
-
+          cl_height = fval(1)
           nz_cl = nz
-
 !
 !--       Find the highest nest level in the coarse grid for the reduced z
@@ -704,5 +687,4 @@
              ENDIF
           ENDDO
-
 !    
 !--       Get absolute coordinates from the child
@@ -711,33 +693,30 @@
          
           CALL pmc_recv_from_child( child_id, cl_coord_x, SIZE( cl_coord_x ),   &
-                                     0, 11, ierr )
+               0, 11, ierr )
           CALL pmc_recv_from_child( child_id, cl_coord_y, SIZE( cl_coord_y ),   &
-                                     0, 12, ierr )
-!          WRITE ( 0, * )  'receive from pmc child ', child_id, nx_cl, ny_cl
+               0, 12, ierr )
          
-          define_coarse_grid_real(1) = lower_left_coord_x
-          define_coarse_grid_real(2) = lower_left_coord_y
-          define_coarse_grid_real(3) = dx
-          define_coarse_grid_real(4) = dy
-          define_coarse_grid_real(5) = lower_left_coord_x + ( nx + 1 ) * dx
-          define_coarse_grid_real(6) = lower_left_coord_y + ( ny + 1 ) * dy
-          define_coarse_grid_real(7) = dz
-
-          define_coarse_grid_int(1)  = nx
-          define_coarse_grid_int(2)  = ny
-          define_coarse_grid_int(3)  = nz_cl
-
+          parent_grid_info_real(1) = lower_left_coord_x
+          parent_grid_info_real(2) = lower_left_coord_y
+          parent_grid_info_real(3) = dx
+          parent_grid_info_real(4) = dy
+          parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx
+          parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy
+          parent_grid_info_real(7) = dz
+
+          parent_grid_info_int(1)  = nx
+          parent_grid_info_int(2)  = ny
+          parent_grid_info_int(3)  = nz_cl
 !
 !--       Check that the child domain matches parent domain. 
           nomatch = 0
           IF ( nesting_mode == 'vertical' )  THEN
-             right_limit = define_coarse_grid_real(5)
-             north_limit = define_coarse_grid_real(6)
+             right_limit = parent_grid_info_real(5)
+             north_limit = parent_grid_info_real(6)
              IF ( ( cl_coord_x(nx_cl+1) /= right_limit ) .OR.                   &
                   ( cl_coord_y(ny_cl+1) /= north_limit ) )  THEN
                 nomatch = 1
              ENDIF
-          ELSE
-          
+          ELSE       
 !
 !--       Check that the child domain is completely inside the parent domain. 
@@ -745,7 +724,7 @@
              yez = ( nbgp + 1 ) * dy 
              left_limit  = lower_left_coord_x + xez
-             right_limit = define_coarse_grid_real(5) - xez
+             right_limit = parent_grid_info_real(5) - xez
              south_limit = lower_left_coord_y + yez
-             north_limit = define_coarse_grid_real(6) - yez
+             north_limit = parent_grid_info_real(6) - yez
              IF ( ( cl_coord_x(0) < left_limit )        .OR.                    &
                   ( cl_coord_x(nx_cl+1) > right_limit ) .OR.                    &
@@ -755,5 +734,11 @@
              ENDIF
           ENDIF
-
+!             
+!--       Child domain must be lower than the parent domain such
+!--       that the top ghost layer of the child grid does not exceed 
+!--       the parent domain top boundary.
+          IF ( cl_height > zw(nz) ) THEN
+             nomatch = 1
+          ENDIF
 !
 !--       Check that parallel nest domains, if any, do not overlap.
@@ -779,13 +764,11 @@
           DEALLOCATE( cl_coord_x )
           DEALLOCATE( cl_coord_y )
-
 !
 !--       Send coarse grid information to child
-          CALL pmc_send_to_child( child_id, define_coarse_grid_real,            &
-                                   SIZE( define_coarse_grid_real ), 0, 21,      &
+          CALL pmc_send_to_child( child_id, parent_grid_info_real,              &
+                                   SIZE( parent_grid_info_real ), 0, 21,        &
                                    ierr )
-          CALL pmc_send_to_child( child_id, define_coarse_grid_int,  3, 0,      &
+          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,        &
                                    22, ierr )
-
 !
 !--       Send local grid to child
@@ -794,5 +777,4 @@
           CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,        &
                                    ierr )
-
 !
 !--       Also send the dzu-, dzw-, zu- and zw-arrays here
@@ -818,9 +800,7 @@
       
        CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
-
 !
 !--    TO_DO: Klaus: please give a comment what is done here
        CALL pmci_create_index_list
-
 !
 !--    Include couple arrays into parent content
@@ -870,4 +850,5 @@
 
        IF ( myid == 0 )  THEN
+!          
 !--       TO_DO: Klaus: give more specific comment what size_of_array stands for
           CALL pmc_recv_from_child( child_id, size_of_array, 2, 0, 40, ierr )
@@ -875,5 +856,4 @@
           CALL pmc_recv_from_child( child_id, coarse_bound_all,                 &
                                     SIZE( coarse_bound_all ), 0, 41, ierr )
-
 !
 !--       Compute size of index_list.
@@ -1023,5 +1003,4 @@
 
        CALL pmc_set_dataarray_name( lastentry = .TRUE. )
-
 !
 !--    Send grid to parent
@@ -1039,37 +1018,34 @@
           CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
           CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
-
 !
 !--       Receive Coarse grid information.
-!--       TO_DO: find shorter and more meaningful name for  define_coarse_grid_real
-          CALL pmc_recv_from_parent( define_coarse_grid_real,                  &
-                                     SIZE(define_coarse_grid_real), 0, 21, ierr )
-          CALL pmc_recv_from_parent( define_coarse_grid_int,  3, 0, 22, ierr )
+          CALL pmc_recv_from_parent( parent_grid_info_real,                     &
+                                     SIZE(parent_grid_info_real), 0, 21, ierr )
+          CALL pmc_recv_from_parent( parent_grid_info_int,  3, 0, 22, ierr )
 !
 !--        Debug-printouts - keep them
 !          WRITE(0,*) 'Coarse grid from parent '
-!          WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
-!          WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
-!          WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(5)
-!          WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(6)
-!          WRITE(0,*) 'dx            = ',define_coarse_grid_real(3)
-!          WRITE(0,*) 'dy            = ',define_coarse_grid_real(4)
-!          WRITE(0,*) 'dz            = ',define_coarse_grid_real(7)
-!          WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
-!          WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
-!          WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)
-       ENDIF
-
-       CALL MPI_BCAST( define_coarse_grid_real, SIZE(define_coarse_grid_real),  &
+!          WRITE(0,*) 'startx_tot    = ',parent_grid_info_real(1)
+!          WRITE(0,*) 'starty_tot    = ',parent_grid_info_real(2)
+!          WRITE(0,*) 'endx_tot      = ',parent_grid_info_real(5)
+!          WRITE(0,*) 'endy_tot      = ',parent_grid_info_real(6)
+!          WRITE(0,*) 'dx            = ',parent_grid_info_real(3)
+!          WRITE(0,*) 'dy            = ',parent_grid_info_real(4)
+!          WRITE(0,*) 'dz            = ',parent_grid_info_real(7)
+!          WRITE(0,*) 'nx_coarse     = ',parent_grid_info_int(1)
+!          WRITE(0,*) 'ny_coarse     = ',parent_grid_info_int(2)
+!          WRITE(0,*) 'nz_coarse     = ',parent_grid_info_int(3)
+       ENDIF
+
+       CALL MPI_BCAST( parent_grid_info_real, SIZE(parent_grid_info_real),      &
                        MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( define_coarse_grid_int, 3, MPI_INTEGER, 0, comm2d, ierr )
-
-       cg%dx = define_coarse_grid_real(3)
-       cg%dy = define_coarse_grid_real(4)
-       cg%dz = define_coarse_grid_real(7)
-       cg%nx = define_coarse_grid_int(1)
-       cg%ny = define_coarse_grid_int(2)
-       cg%nz = define_coarse_grid_int(3)
-
+       CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr )
+
+       cg%dx = parent_grid_info_real(3)
+       cg%dy = parent_grid_info_real(4)
+       cg%dz = parent_grid_info_real(7)
+       cg%nx = parent_grid_info_int(1)
+       cg%ny = parent_grid_info_int(2)
+       cg%nz = parent_grid_info_int(3)
 !
 !--    Get parent coordinates on coarse grid
@@ -1081,5 +1057,4 @@
        ALLOCATE( cg%zu(0:cg%nz+1) )
        ALLOCATE( cg%zw(0:cg%nz+1) )
-
 !
 !--    Get coarse grid coordinates and values of the z-direction from the parent
@@ -1094,5 +1069,4 @@
 
        ENDIF
-
 !
 !--    Broadcast this information
@@ -1110,5 +1084,4 @@
 !--    TO_DO: Klaus give a comment what is happening here
        CALL pmc_c_get_2d_index_list
-
 !
 !--    Include couple arrays into child content
@@ -1116,21 +1089,18 @@
        CALL  pmc_c_clear_next_array_list
        DO  WHILE ( pmc_c_getnextarray( myname ) )
-!--       TO_DO: Klaus, why the child-arrays are still up to cg%nz??
+!--       Note that cg%nz is not th eoriginal nz of parent, but the highest
+!--       parent-grid level needed for nesting.           
           CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz )
        ENDDO
        CALL pmc_c_setind_and_allocmem
-
 !
 !--    Precompute interpolation coefficients and child-array indices
        CALL pmci_init_interp_tril
-
 !
 !--    Precompute the log-law correction index- and ratio-arrays 
        CALL pmci_init_loglaw_correction 
-
 !
 !--    Define the SGS-TKE scaling factor based on the grid-spacing ratio
        CALL pmci_init_tkefactor
-
 !
 !--    Two-way coupling for general and vertical nesting.
@@ -1141,5 +1111,4 @@
           CALL pmci_init_anterp_tophat
        ENDIF
-
 !
 !--    Finally, compute the total area of the top-boundary face of the domain.
@@ -1150,4 +1119,5 @@
 
  CONTAINS
+
 
     SUBROUTINE pmci_map_fine_to_coarse_grid
@@ -1263,4 +1233,6 @@
        INTEGER(iwp) ::  k       !:
        INTEGER(iwp) ::  kc      !:
+       INTEGER(iwp) ::  kdzo    !:
+       INTEGER(iwp) ::  kdzw    !:       
 
        REAL(wp) ::  xb          !:
@@ -1301,5 +1273,4 @@
        ALLOCATE( r1zo(nzb:nzt+1) )
        ALLOCATE( r2zo(nzb:nzt+1) )
-
 !
 !--    Note that the node coordinates xfs... and xcs... are relative to the
@@ -1336,13 +1307,18 @@
           zfso = zu(k)
 
-          kc = 0
-          DO  WHILE ( cg%zw(kc) <= zfsw )
-             kc = kc + 1
+          DO kc = 0, cg%nz+1
+             IF ( cg%zw(kc) > zfsw )  EXIT
           ENDDO
           kcw(k) = kc - 1
-         
-          kc = 0
-          DO  WHILE ( cg%zu(kc) <= zfso )
-             kc = kc + 1
+          
+          !if ( myid == 0 .and. nx==191 )  then
+          !   write(162,*)nx, nzt+1, k, zw(k), cg%nz+1, kcw(k), cg%zw(kcw(k))
+          !endif
+          !if ( myid == 0 .and. nx==383 )  then
+          !   write(163,*)nx, nzt+1, k, zw(k), cg%nz+1, kcw(k), cg%zw(kcw(k))
+          !endif
+
+          DO kc = 0, cg%nz+1
+             IF ( cg%zu(kc) > zfso )  EXIT
           ENDDO
           kco(k) = kc - 1
@@ -1350,6 +1326,8 @@
           zcsw    = cg%zw(kcw(k))
           zcso    = cg%zu(kco(k))
-          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k)+1)
-          r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k)+1)
+          kdzw    = MIN( kcw(k)+1, cg%nz+1 )
+          kdzo    = MIN( kco(k)+1, cg%nz+1 )
+          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kdzw)
+          r2zo(k) = ( zfso - zcso ) / cg%dzu(kdzo)
           r1zw(k) = 1.0_wp - r2zw(k)
           r1zo(k) = 1.0_wp - r2zo(k)
@@ -1507,5 +1485,4 @@
           nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
        ENDIF
-
 !
 !--    Then determine the maximum number of near-wall nodes per wall point based
@@ -1513,5 +1490,4 @@
        nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,                &
                            nzt_topo_nestbc_s, nzt_topo_nestbc_n )
-
 !
 !--    Note that the outer division must be integer division.
@@ -1529,5 +1505,4 @@
 
        z0_topo = roughness_length
-
 !
 !--    First horizontal walls. Note that also logc_w_? and logc_ratio_w_? need to 
@@ -1590,5 +1565,4 @@
 
        ENDIF
-
 !
 !--    Right boundary
@@ -1650,5 +1624,4 @@
 
        ENDIF
-
 !
 !--    South boundary
@@ -1707,5 +1680,4 @@
 
        ENDIF
-
 !
 !--    North boundary
@@ -1763,5 +1735,4 @@
 
        ENDIF
-
 !       
 !--    Then vertical walls and corners if necessary
@@ -1817,5 +1788,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for u on the north side, but not on the south side
@@ -1851,5 +1821,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the north side, but not on the south side.
@@ -1872,5 +1841,4 @@
 
           ENDIF   !  IF ( nest_bound_l )
-
 !        
 !--       Right boundary
@@ -1907,5 +1875,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for u on the north side, but not on the south side
@@ -1923,5 +1890,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the south side, but not on the north side
@@ -1939,5 +1905,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the north side, but not on the south side
@@ -1948,5 +1913,4 @@
                       CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
                           k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
-
 !
 !--                   The direction of the wall-normal index is stored as the
@@ -1961,5 +1925,4 @@
 
           ENDIF   !  IF ( nest_bound_r )
-
 !        
 !--       South boundary
@@ -1996,5 +1959,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for v on the right side, but not on the left side
@@ -2013,5 +1975,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the left side, but not on the right side
@@ -2052,5 +2013,4 @@
 
           ENDIF   !  IF (nest_bound_s )
-
 !        
 !--       North boundary
@@ -2087,5 +2047,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for v on the right side, but not on the left side
@@ -2103,5 +2062,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the left side, but not on the right side
@@ -2119,5 +2077,4 @@
                       lcr(0:ncorr-1) = 1.0_wp
                    ENDIF
-
 !
 !--                Wall for w on the right side, but not on the left side
@@ -2204,5 +2161,4 @@
                                                 z0_l, inc )
                 ENDIF
-
 !
 !--             The role of inc here is to make the comparison operation "<"
@@ -2311,5 +2267,5 @@
        REAL(wp) ::  logyc1   !:
        REAL(wp) ::  yc1      !:
-
+       
 !
 !--    yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from
@@ -2416,5 +2372,4 @@
           anterp_relax_length_t = 0.1_wp * zu(nzt)
        ENDIF
-
 !
 !--    First determine kctu and kctw that are the coarse-grid upper bounds for
@@ -2472,5 +2427,4 @@
        iflu(icr) = nxrg
        ifuu(icr) = nxrg
-
 !
 !--    i-indices of others for each ii-index value
@@ -2495,5 +2449,4 @@
        iflo(icr) = nxrg
        ifuo(icr) = nxrg
-
 !
 !--    j-indices of v for each jj-index value
@@ -2540,5 +2493,4 @@
        jflo(jcn) = nyng
        jfuo(jcn) = nyng
-
 !
 !--    k-indices of w for each kk-index value
@@ -2558,5 +2510,4 @@
           kstart = kflw(kk)
        ENDDO
-
 !
 !--    k-indices of others for each kk-index value
@@ -2576,5 +2527,4 @@
           kstart = kflo(kk)
        ENDDO
- 
 !
 !--    Precomputation of number of fine-grid nodes inside coarse-grid ij-faces.
@@ -2645,5 +2595,5 @@
 
 
-SUBROUTINE pmci_init_tkefactor
+    SUBROUTINE pmci_init_tkefactor
 
 !
@@ -2712,5 +2662,5 @@
        ENDIF
 
-      IF ( nest_bound_s )  THEN
+       IF ( nest_bound_s )  THEN
           ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
           tkefactor_s = 0.0_wp
@@ -2718,7 +2668,7 @@
           DO  i = nxlg, nxrg
              k_wall = get_topography_top_index( j, i, 's' )
-
+             
              DO  k = k_wall + 1, nzt
-
+                
                 kc     = kco(k+1)
                 glsf   = ( dx * dy * dzu(k) )**p13
@@ -2727,5 +2677,5 @@
                 fw     = EXP( -cfw*height / glsf )
                 tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
-                                              ( glsf / glsc )**p23 )
+                     ( glsf / glsc )**p23 )
              ENDDO
              tkefactor_s(k_wall,i) = c_tkef * fw0
@@ -2827,5 +2777,4 @@
     NULLIFY( p_3d )
     NULLIFY( p_2d )
-
 !
 !-- List of array names, which can be coupled.
@@ -2928,5 +2877,4 @@
     NULLIFY( p_3d )
     NULLIFY( p_2d )
-
 !
 !-- List of array names, which can be coupled
@@ -3036,7 +2984,6 @@
 
 !
-!-- Root id is never a child
+!-- Root model is never anyone's child
     IF ( cpl_id > 1 )  THEN
-
 !
 !--    Child domain boundaries in the parent index space
@@ -3045,9 +2992,7 @@
        jcs = coarse_bound(3)
        jcn = coarse_bound(4)
-
 !
 !--    Get data from the parent
        CALL pmc_c_getbuffer( waittime = waittime )
-
 !
 !--    The interpolation.
@@ -3203,5 +3148,4 @@
           ENDIF
        ENDIF
-
 !
 !--    Trilinear interpolation.
@@ -3222,5 +3166,4 @@
           ENDDO
        ENDDO
-
 !
 !--    Correct the interpolated values of u and v in near-wall nodes, i.e. in
@@ -3321,5 +3264,4 @@
        ENDIF
     ENDIF
-
 !
 !-- Same for restart time
@@ -3338,5 +3280,4 @@
        ENDIF
     ENDIF
-
 !
 !-- Same for dt_restart
@@ -3355,5 +3296,4 @@
        ENDIF
     ENDIF
-
 !
 !-- Same for time_restart
@@ -3432,5 +3372,5 @@
     volume_flow(1) = volume_flow_l(1)
 #endif
-
+    
 !
 !-- Sum up the volume flow through the south/north boundaries
@@ -3512,4 +3452,5 @@
 !-- MPI_ALLREDUCE with the MPI_MIN operator over all processes using 
 !-- the global communicator MPI_COMM_WORLD.
+   
    IMPLICIT NONE
 
@@ -3605,5 +3546,4 @@
 
  END SUBROUTINE pmci_datatrans
-
 
 
@@ -3642,9 +3582,7 @@
           CALL pmc_s_getdata_from_buffer( child_id )
           CALL cpu_log( log_point_s(72), 'pmc parent recv', 'stop' )
-
 !
 !--       The anterpolated data is now available in u etc
           IF ( topography /= 'flat' )  THEN
-
 !
 !--          Inside buildings/topography reset velocities back to zero.
@@ -3732,12 +3670,14 @@
     ENDIF
 
- CONTAINS
-
+  CONTAINS
+
+   
     SUBROUTINE pmci_interpolation
 
 !
 !--    A wrapper routine for all interpolation and extrapolation actions
+      
        IMPLICIT NONE
-
+      
 !
 !--    In case of vertical nesting no interpolation is needed for the
@@ -3746,5 +3686,5 @@
        
 !
-!--    Left border pe:
+!--       Left border pe:
           IF ( nest_bound_l )  THEN
              CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
@@ -3855,7 +3795,6 @@
 
           ENDIF
-
-   !
-   !--    Right border pe
+!
+!--       Right border pe
           IF ( nest_bound_r )  THEN
              CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
@@ -3969,7 +3908,6 @@
 
           ENDIF
-
-   !
-   !--    South border pe
+!
+!--       South border pe
           IF ( nest_bound_s )  THEN
              CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
@@ -4077,7 +4015,6 @@
 
           ENDIF
-
-   !
-   !--    North border pe
+!
+!--       North border pe
           IF ( nest_bound_n )  THEN
              CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
@@ -4189,5 +4126,4 @@
 
        ENDIF       ! IF ( nesting_mode /= 'vertical' )
-
 !
 !--    All PEs are top-border PEs
@@ -4413,5 +4349,9 @@
       
       DO  j = nys, nyn+1
-         DO  k = nzb, nzt+1
+!
+!--      Determine vertical index of topography top at grid point (j,i)
+         k_wall = get_topography_top_index( j, i, TRIM( var ) )
+
+         DO  k = k_wall, nzt+1
             l = ic(i)
             m = jc(j)
@@ -4426,5 +4366,4 @@
          ENDDO
       ENDDO
-
 !
 !--   Generalized log-law-correction algorithm.
@@ -4450,5 +4389,4 @@
          ENDDO
       ENDIF
-
 !
 !--   In case of non-flat topography, also vertical walls and corners need to be
@@ -4513,5 +4451,4 @@
 
       ENDIF  ! ( topography /= 'flat' )
-
 !
 !--   Rescale if f is the TKE.
@@ -4539,5 +4476,4 @@
          ENDIF
       ENDIF
-
 !
 !--   Store the boundary values also into the other redundant ghost node layers
@@ -4634,5 +4570,4 @@
       ENDIF
 
-
       DO  i = nxl, nxr+1
 !
@@ -4653,5 +4588,4 @@
          ENDDO
       ENDDO
-
 !
 !--   Generalized log-law-correction algorithm.
@@ -4677,5 +4611,4 @@
          ENDDO
       ENDIF
-
 !
 !--   In case of non-flat topography, also vertical walls and corners need to be
@@ -4741,5 +4674,4 @@
          
       ENDIF  ! ( topography /= 'flat' )
-
 !
 !--   Rescale if f is the TKE.
@@ -4765,5 +4697,4 @@
          ENDIF
       ENDIF
-
 !
 !--   Store the boundary values also into the other redundant ghost node layers
@@ -4810,5 +4741,9 @@
 
       INTEGER(iwp) ::  i   !:
+      INTEGER(iwp) ::  ib  !:
+      INTEGER(iwp) ::  ie  !:
       INTEGER(iwp) ::  j   !:
+      INTEGER(iwp) ::  jb   !:
+      INTEGER(iwp) ::  je   !:      
       INTEGER(iwp) ::  k   !:
       INTEGER(iwp) ::  l   !:
@@ -4832,10 +4767,27 @@
          k  = nzt + 1
       ENDIF
-     
-      DO  i = nxl-1, nxr+1
-         DO  j = nys-1, nyn+1
+!
+!--   These exceedings by one are needed only to avoid stripes
+!--   and spots in visualization. They have no effect on the 
+!--   actual solution.     
+      ib = nxl-1
+      ie = nxr+1
+      jb = nys-1
+      je = nyn+1
+!
+!--   The exceedings must not be made past the outer edges in 
+!--   case of pure vertical nesting.
+      IF ( nesting_mode == 'vertical' )  THEN
+         IF ( nxl == 0  )  ib = nxl
+         IF ( nxr == nx )  ie = nxr
+         IF ( nys == 0  )  jb = nys
+         IF ( nyn == ny )  je = nyn
+      ENDIF
+         
+      DO  i = ib, ie
+         DO  j = jb, je
             l = ic(i)
             m = jc(j)
-            n = kc(k)              
+            n = kc(k)            
             fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
             fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
@@ -4847,5 +4799,4 @@
          ENDDO
       ENDDO
-
 !
 !--   Just fill up the second ghost-node layer for w.
@@ -4853,5 +4804,4 @@
          f(nzt+1,:,:) = f(nzt,:,:)
       ENDIF
-
 !
 !--   Rescale if f is the TKE.
@@ -4896,5 +4846,4 @@
 
        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
-
 !
 !--    Check which edge is to be handled: left or right
@@ -4917,5 +4866,4 @@
        ENDIF
 
-
        DO  j = nys, nyn+1
 !
@@ -4935,5 +4883,4 @@
           ENDIF
        ENDDO
-
 !
 !--    Store the boundary values also into the redundant ghost node layers.
@@ -4998,5 +4945,4 @@
        ENDIF
 
-
        DO  i = nxl, nxr+1
 !
@@ -5016,5 +4962,4 @@
           ENDIF
        ENDDO
-
 !
 !--    Store the boundary values also into the redundant ghost node layers.
@@ -5072,5 +5017,4 @@
           ENDDO
        ENDDO
-
 !
 !--    Just fill up the second ghost-node layer for w
@@ -5126,5 +5070,4 @@
        REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)  ::  fc  !< Treated variable - parent domain
   
-
 !
 !--    Initialize the index bounds for anterpolation 
@@ -5179,6 +5122,5 @@
              jcnm = jcn - nhln - 1
           ENDIF
-       ENDIF
-       
+       ENDIF     
 !
 !--    Note that ii, jj, and kk are coarse-grid indices and i,j, and k 
@@ -5211,5 +5153,4 @@
           ENDDO
        ENDDO
-
 
     END SUBROUTINE pmci_anterp_tophat
@@ -5298,5 +5239,4 @@
        ENDIF
     ENDIF
-
 !
 !-- Set Neumann boundary conditions for humidity and cloud-physical quantities
