Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 3483)
+++ /palm/trunk/SOURCE/palm.f90	(revision 3484)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! pmci_ensure_nest_mass_conservation removed permanently
+! 
+! 3458 2018-10-30 14:51:23Z kanani
 ! from chemistry branch r3443, forkel:
 ! removed double do_emis check around CALL chem_init
@@ -308,6 +311,5 @@
     USE pmc_interface,                                                         &
         ONLY:  nested_run, pmci_child_initialize, pmci_init,                   &
-               pmci_modelconfiguration, pmci_parent_initialize,                &
-               pmci_ensure_nest_mass_conservation
+               pmci_modelconfiguration, pmci_parent_initialize
 
     USE write_restart_data_mod,                                                &
Index: /palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3483)
+++ /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3484)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Introduction of reversibility correction to the interpolation routines in order to 
+! guarantee mass and scalar conservation through the nest boundaries. Several errors 
+! are corrected and pmci_ensure_nest_mass_conservation is permanently removed. 
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -366,5 +371,5 @@
     USE pegrid,                                                                &
         ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
-               numprocs
+               numprocs, pleft, pnorth, pright, psouth, status
 
     USE pmc_child,                                                             &
@@ -417,5 +422,5 @@
     INTEGER(iwp), SAVE      ::  cpl_parent_id         !<
 !
-!-- Control parameters, will be made input parameters later
+!-- Control parameters
     CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !< steering
                                                          !< parameter for data-
@@ -441,6 +446,7 @@
 
 !
-!-- Child coarse data arrays
-    INTEGER(iwp), DIMENSION(5),PUBLIC           ::  coarse_bound   !<
+!-- Children's parent-grid arrays
+    INTEGER(iwp), SAVE, DIMENSION(5), PUBLIC    ::  coarse_bound        !< subdomain index bounds for children's parent-grid arrays
+    INTEGER(iwp), SAVE, DIMENSION(4), PUBLIC    ::  coarse_bound_anterp !< subdomain index bounds for anterpolation
 
     REAL(wp), SAVE                              ::  xexl           !<
@@ -492,4 +498,5 @@
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !<
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !<
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  celltmpd !<
 !
 !-- Child index arrays and log-ratio arrays for the log-law near-wall
@@ -534,5 +541,5 @@
 !
 !-- Upper bounds for k in anterpolation.
-    INTEGER(iwp), SAVE ::  kctu   !<
+    INTEGER(iwp), SAVE ::  kcto   !<
     INTEGER(iwp), SAVE ::  kctw   !<
 !
@@ -599,7 +606,6 @@
 
     TYPE(coarsegrid_def), SAVE, PUBLIC     ::  cg   !<
-
-!-  Variables for particle coupling
-
+!
+!-- Variables for particle coupling
     TYPE, PUBLIC :: childgrid_def
        INTEGER(iwp)                        ::  nx                   !<
@@ -641,8 +647,4 @@
     END INTERFACE pmci_datatrans
 
-    INTERFACE pmci_ensure_nest_mass_conservation
-       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
-    END INTERFACE
-
     INTERFACE pmci_init
        MODULE PROCEDURE pmci_init
@@ -687,5 +689,4 @@
     PUBLIC pmci_child_initialize
     PUBLIC pmci_datatrans
-    PUBLIC pmci_ensure_nest_mass_conservation
     PUBLIC pmci_init
     PUBLIC pmci_modelconfiguration
@@ -891,6 +892,6 @@
        IF ( myid == 0 )  THEN
 
-          CALL pmc_recv_from_child( child_id, val,  size(val),  0, 123, ierr )
-          CALL pmc_recv_from_child( child_id, fval, size(fval), 0, 124, ierr )
+          CALL pmc_recv_from_child( child_id, val,  SIZE(val),  0, 123, ierr )
+          CALL pmc_recv_from_child( child_id, fval, SIZE(fval), 0, 124, ierr )
          
           nx_cl     = val(1)
@@ -983,5 +984,5 @@
                    mid = pmc_parent_for_child(mm)
 !
-!--                Check Only different nest level
+!--                Check only different nest levels
                    IF (m_couplers(child_id)%parent_id /= m_couplers(mid)%parent_id)  THEN
                       IF ( ( ch_xl(m) < ch_xr(mm) .OR.                         &
@@ -1000,5 +1001,4 @@
           DEALLOCATE( cl_coord_x )
           DEALLOCATE( cl_coord_y )
-
 !
 !--       Send information about operating mode (LES or RANS) to child. This will be 
@@ -1219,10 +1219,14 @@
 
     INTEGER(iwp) ::  i          !<
-    INTEGER(iwp) ::  ierr       !<
-    INTEGER(iwp) ::  icl        !<
-    INTEGER(iwp) ::  icr        !<
-    INTEGER(iwp) ::  j          !<
-    INTEGER(iwp) ::  jcn        !<
-    INTEGER(iwp) ::  jcs        !<
+    INTEGER(iwp) ::  ierr       !< 
+    INTEGER(iwp) ::  icl        !< left index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  icla       !< left index limit for anterpolation
+    INTEGER(iwp) ::  icr        !< left index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  icra       !< right index limit for anterpolation
+    INTEGER(iwp) ::  j          !< 
+    INTEGER(iwp) ::  jcn        !< north index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  jcna       !< north index limit for anterpolation
+    INTEGER(iwp) ::  jcs        !< sout index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  jcsa       !< south index limit for anterpolation
     INTEGER(iwp) ::  n          !< running index for number of chemical species
 
@@ -1374,12 +1378,10 @@
 !--    Get coarse grid coordinates and values of the z-direction from the parent
        IF ( myid == 0)  THEN
-
           CALL pmc_recv_from_parent( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr )
           CALL pmc_recv_from_parent( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr )
-          CALL pmc_recv_from_parent( cg%dzu, cg%nz + 1, 0, 26, ierr )
-          CALL pmc_recv_from_parent( cg%dzw, cg%nz + 1, 0, 27, ierr )
-          CALL pmc_recv_from_parent( cg%zu, cg%nz + 2, 0, 28, ierr )
-          CALL pmc_recv_from_parent( cg%zw, cg%nz + 2, 0, 29, ierr )
-
+          CALL pmc_recv_from_parent( cg%dzu, cg%nz+1, 0, 26, ierr )
+          CALL pmc_recv_from_parent( cg%dzw, cg%nz+1, 0, 27, ierr )
+          CALL pmc_recv_from_parent( cg%zu, cg%nz+2, 0, 28, ierr )
+          CALL pmc_recv_from_parent( cg%zw, cg%nz+2, 0, 29, ierr )
        ENDIF
 !
@@ -1392,5 +1394,4 @@
        CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
        CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )
-  
 !
 !--    Find the index bounds for the nest domain in the coarse-grid index space
@@ -1441,8 +1442,9 @@
 !--    Precompute the index arrays and relaxation functions for the
 !--    anterpolation
-       IF ( TRIM( nesting_mode ) == 'two-way' .OR.                             &
-                  nesting_mode == 'vertical' )  THEN
-          CALL pmci_init_anterp_tophat
-       ENDIF
+!
+!--    Note that the anterpolation index bounds are needed also in case
+!--    of one-way coupling because of the reversibility correction
+!--    included in the interpolation algorithms.
+       CALL pmci_init_anterp_tophat
 !
 !--    Finally, compute the total area of the top-boundary face of the domain.
@@ -1463,5 +1465,8 @@
        INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all   !<
        INTEGER(iwp), DIMENSION(2)          ::  size_of_array      !<
-                                              
+                                   
+       INTEGER(iwp) :: i        !<   
+       INTEGER(iwp) :: ijaux    !< 
+       INTEGER(iwp) :: j        !< 
        REAL(wp) ::  loffset     !<
        REAL(wp) ::  noffset     !<
@@ -1536,4 +1541,64 @@
        coarse_bound(5) = myid
 !
+!--    Determine the anterpolation index limits. If at least half of the 
+!--    parent-grid cell is within the current child sub-domain, then it 
+!--    is included in the current sub-domain's anterpolation domain. 
+!--    Else the parent-grid cell is included in the neighbouring subdomain's 
+!--    anterpolation domain, or not included at all if we are at the outer 
+!--    edge of the child domain.
+       DO  i = 0, cg%nx
+          IF ( cg%coord_x(i) + 0.5_wp * cg%dx >= coord_x(nxl) )  THEN
+             icla = MAX( 0, i )
+             EXIT
+          ENDIF
+       ENDDO
+       DO  i = cg%nx, 0 , -1
+          IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= coord_x(nxr+1) )  THEN
+             icra = MIN( cg%nx, i )
+             EXIT
+          ENDIF
+       ENDDO
+       DO  j = 0, cg%ny
+          IF ( cg%coord_y(j) + 0.5_wp * cg%dy >= coord_y(nys) )  THEN
+             jcsa = MAX( 0, j )
+             EXIT
+          ENDIF
+       ENDDO
+       DO  j = cg%ny, 0 , -1
+          IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= coord_y(nyn+1) )  THEN
+             jcna = MIN( cg%ny, j )
+             EXIT
+          ENDIF
+       ENDDO
+!
+!--    Make sure that the indexing is contiguous 
+       IF ( nxl == 0 )  THEN
+          CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+       ELSE IF ( nxr == nx )  THEN
+          CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr )
+          icla = ijaux + 1
+       ELSE
+          CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+          CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 
+          icla = ijaux + 1
+       ENDIF
+       IF ( nys == 0 )  THEN
+          CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+       ELSE IF ( nyn == ny )  THEN
+          CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr )
+          jcsa = ijaux + 1
+       ELSE
+          CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+          CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 
+          jcsa = ijaux + 1
+       ENDIF
+
+       write(9,"('Anterpolation bounds: ',4(i3,2x))") icla, icra, jcsa, jcna
+       flush(9)
+       coarse_bound_anterp(1) = icla
+       coarse_bound_anterp(2) = icra
+       coarse_bound_anterp(3) = jcsa
+       coarse_bound_anterp(4) = jcna
+!
 !--    Note that MPI_Gather receives data from all processes in the rank order
 !--    TO_DO: refer to the line where this fact becomes important
@@ -1561,4 +1626,5 @@
        IMPLICIT NONE
 
+       INTEGER(iwp) ::  acsize  !<
        INTEGER(iwp) ::  i       !<
        INTEGER(iwp) ::  j       !<
@@ -1568,4 +1634,6 @@
        INTEGER(iwp) ::  kdzw    !<       
 
+       REAL(wp) ::  dzmin       !<
+       REAL(wp) ::  parentdzmax !<
        REAL(wp) ::  xb          !<
        REAL(wp) ::  xcsu        !<
@@ -1583,5 +1651,5 @@
        REAL(wp) ::  zfsw        !<
       
-
+       
        xb = nxl * dx
        yb = nys * dy
@@ -1658,4 +1726,20 @@
           r1zo(k) = 1.0_wp - r2zo(k)
        ENDDO
+!
+!--    Determine the maximum dimension of anterpolation cells and allocate the
+!--    work array celltmpd needed in the reversibility correction in the 
+!--    interpolation
+       dzmin = 999999.9_wp
+       DO k = 1, nzt+1
+          dzmin = MIN( dzmin, dzu(k), dzw(k) )
+       ENDDO
+       parentdzmax = 0.0_wp
+       DO k = 1, cg%nz+1
+          parentdzmax = MAX(parentdzmax , cg%dzu(k), cg%dzw(k) )
+       ENDDO
+       acsize = CEILING( cg%dx / dx ) * CEILING( cg%dy / dy ) *                 &
+            CEILING( parentdzmax / dzmin )
+       ALLOCATE( celltmpd(1:acsize) )
+!       write(9,"('acsize: ',i3,2(e12.5,2x))") acsize, dzmin, parentdzmax
       
     END SUBROUTINE pmci_init_interp_tril
@@ -3007,4 +3091,5 @@
        REAL(wp)     ::  xi       !<
        REAL(wp)     ::  eta      !<
+       REAL(wp)     ::  tolerance !<
        REAL(wp)     ::  zeta     !<
      
@@ -3027,5 +3112,5 @@
        ENDIF
 !
-!--    First determine kctu and kctw that are the coarse-grid upper bounds for
+!--    First determine kcto and kctw that are the coarse-grid upper bounds for
 !--    index k
        kk = 0
@@ -3033,5 +3118,5 @@
           kk = kk + 1
        ENDDO
-       kctu = kk - 1
+       kcto = kk - 1
 
        kk = 0
@@ -3049,13 +3134,14 @@
        ALLOCATE( jfuv(jcs:jcn) )
        ALLOCATE( jfuo(jcs:jcn) )
-       ALLOCATE( kflw(0:kctw) )
-       ALLOCATE( kflo(0:kctu) )
-       ALLOCATE( kfuw(0:kctw) )
-       ALLOCATE( kfuo(0:kctu) )
-
-       ALLOCATE( ijkfc_u(0:kctu,jcs:jcn,icl:icr) )
-       ALLOCATE( ijkfc_v(0:kctu,jcs:jcn,icl:icr) )
-       ALLOCATE( ijkfc_w(0:kctw,jcs:jcn,icl:icr) )
-       ALLOCATE( ijkfc_s(0:kctu,jcs:jcn,icl:icr) )
+       ALLOCATE( kflw(0:cg%nz+1) )
+       ALLOCATE( kflo(0:cg%nz+1) )
+       ALLOCATE( kfuw(0:cg%nz+1) )
+       ALLOCATE( kfuo(0:cg%nz+1) )
+
+       ALLOCATE( ijkfc_u(0:cg%nz+1,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_v(0:cg%nz+1,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_w(0:cg%nz+1,jcs:jcn,icl:icr) )
+       ALLOCATE( ijkfc_s(0:cg%nz+1,jcs:jcn,icl:icr) )
+
        ijkfc_u = 0
        ijkfc_v = 0
@@ -3065,20 +3151,44 @@
 !--    i-indices of u for each ii-index value
 !--    ii=icr is redundant for anterpolation
+       tolerance = 0.000001_wp * dx
        istart = nxlg
        DO  ii = icl, icr-1
+!
+!--       In case the child and parent grid lines match in x 
+!--       use only the local k,j-child-grid plane for the anterpolation, 
+!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
           i = istart
-          DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  &
-                      ( i < nxrg ) )
-             i  = i + 1
+          DO WHILE ( coord_x(i) < cg%coord_x(ii) - tolerance  .AND.  i < nxrg )
+             i = i + 1
           ENDDO
-          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
-          ir = i
-          DO  WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.&
-                      ( i < nxrg+1 ) )
-             i  = i + 1
-             ir = MIN( i, nxrg )
-          ENDDO
-          ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg )
-          istart = iflu(ii)
+          IF ( ABS( coord_x(i) - cg%coord_x(ii) ) < tolerance )  THEN
+             i = istart
+             DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg )
+                i = i + 1
+             ENDDO
+             iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
+             ifuu(ii) = iflu(ii)
+             istart   = iflu(ii)
+          ELSE
+             i = istart
+             DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )       &
+                  .AND.  ( i < nxrg ) )
+                i  = i + 1
+             ENDDO
+             iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
+             ir = i
+             DO  WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx )     &
+                  .AND.  ( i < nxrg+1 ) )
+                i  = i + 1
+                ir = MIN( i, nxrg )
+             ENDDO
+             ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg )
+             istart = iflu(ii)
+          ENDIF
+!AH
+          write(9,"('pmci_init_anterp_tophat, ii, iflu, ifuu: ', 3(i4,2x))")    &
+               ii, iflu(ii), ifuu(ii)
+          flush(9)
+
        ENDDO
        iflu(icr) = nxrg
@@ -3104,4 +3214,9 @@
           istart = iflo(ii)
        ENDDO
+!AH
+       write(9,"('pmci_init_anterp_tophat, ii, iflo, ifuo: ', 3(i4,2x))")    &
+            ii, iflo(ii), ifuo(ii)
+       flush(9)
+          
        iflo(icr) = nxrg
        ifuo(icr) = nxrg
@@ -3109,20 +3224,44 @@
 !--    j-indices of v for each jj-index value
 !--    jj=jcn is redundant for anterpolation
+       tolerance = 0.000001_wp * dy
        jstart = nysg
        DO  jj = jcs, jcn-1
+!
+!--       In case the child and parent grid lines match in y 
+!--       use only the local k,i-child-grid plane for the anterpolation, 
+!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
           j = jstart
-          DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  &
-                      ( j < nyng ) )
-             j  = j + 1
+          DO WHILE ( coord_y(j) < cg%coord_y(jj) - tolerance  .AND.  j < nyng )
+             j = j + 1
           ENDDO
-          jflv(jj) = MIN( MAX( j, nysg ), nyng )
-          jr = j
-          DO  WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.&
-                      ( j < nyng+1 ) )
-             j  = j + 1
-             jr = MIN( j, nyng )
-          ENDDO
-          jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng )
-          jstart = jflv(jj)
+          IF ( ABS( coord_y(j) - cg%coord_y(jj) ) < tolerance )  THEN
+             j = jstart
+             DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng )
+                j = j + 1
+             ENDDO
+             jflv(jj) = MIN( MAX( j, nysg ), nyng )
+             jfuv(jj) = jflv(jj)
+             jstart   = jflv(jj)
+          ELSE
+             j = jstart
+             DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )       &
+                  .AND.  ( j < nyng ) )
+                j  = j + 1
+             ENDDO
+             jflv(jj) = MIN( MAX( j, nysg ), nyng )
+             jr = j
+             DO  WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy )     &
+                  .AND.  ( j < nyng+1 ) )
+                j  = j + 1
+                jr = MIN( j, nyng )
+             ENDDO
+             jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng )
+             jstart = jflv(jj)
+          ENDIF
+!AH
+          write(9,"('pmci_init_anterp_tophat, jj, jflv, jfuv: ', 3(i4,2x))")    &
+               jj, jflv(jj), jfuv(jj)
+          flush(9)
+
        ENDDO
        jflv(jcn) = nyng
@@ -3147,4 +3286,9 @@
           jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng )
           jstart = jflo(jj)
+!AH
+          write(9,"('pmci_init_anterp_tophat, ii, jflo, jfuo: ', 3(i4,2x))")    &
+               jj, jflo(jj), jfuo(jj)
+          flush(9)
+
        ENDDO
        jflo(jcn) = nyng
@@ -3152,18 +3296,49 @@
 !
 !--    k-indices of w for each kk-index value
+!--    Note that anterpolation index limits are needed also for the top boundary 
+!--    ghost cell level because of the reversibility correction in the interpolation.
        kstart  = 0
        kflw(0) = 0
        kfuw(0) = 0
-       DO  kk = 1, kctw
+       tolerance = 0.000001_wp * dzw(1)
+       DO kk = 1, cg%nz+1
+!
+!--       In case the child and parent grid lines match in z
+!--       use only the local j,i-child-grid plane for the anterpolation, 
+!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
           k = kstart
-          DO  WHILE ( ( zw(k) < cg%zu(kk) )  .AND.  ( k < nzt ) )
+          DO WHILE ( zw(k) < cg%zw(kk) - tolerance  .AND.  k < nzt+1 )
              k = k + 1
           ENDDO
-          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-          DO  WHILE ( ( zw(k) <= cg%zu(kk+1) )  .AND.  ( k < nzt+1 ) )
-             k  = k + 1
-          ENDDO
-          kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 )
-          kstart = kflw(kk)
+          IF ( ABS( zw(k) - cg%zw(kk) ) < tolerance )  THEN
+             k = kstart
+             DO WHILE ( ( zw(k) < cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
+                k = k + 1
+             ENDDO
+             kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+             kfuw(kk) = kflw(kk)
+             kstart   = kflw(kk)
+          ELSE
+             k = kstart
+             DO  WHILE ( ( zw(k) < cg%zu(kk) )  .AND.  ( k < nzt ) )
+                k = k + 1
+             ENDDO
+             kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+             IF ( kk+1 <= cg%nz+1 )  THEN
+                DO  WHILE ( ( zw(k) <= cg%zu(kk+1) )  .AND.  ( k < nzt+1 ) )
+                   k  = k + 1
+                   IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
+                ENDDO
+                kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 )
+             ELSE
+                kfuw(kk) = kflw(kk)
+             ENDIF
+             kstart = kflw(kk)
+          ENDIF
+!AH
+          write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 3(i4,2x))") &
+               kk, kflw(kk), kfuw(kk)
+          flush(9)
+
        ENDDO
 !
@@ -3172,23 +3347,36 @@
        kflo(0) = 0
        kfuo(0) = 0
-       DO  kk = 1, kctu
+!
+!--    Note that anterpolation index limits are needed also for the top boundary 
+!--    ghost cell level because of the reversibility correction in the interpolation.
+!AH       DO  kk = 1, kcto+1
+       DO  kk = 1, cg%nz+1
           k = kstart
-          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k < nzt ) )
+!AH          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k < nzt ) )
+!--       Note that this is an IMPORTANT correction for the reversibility correction
+          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k <= nzt ) )
              k = k + 1
           ENDDO
           kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
+!AH          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
+!--       Note that this is an IMPORTANT correction for the reversibility correction
+          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
              k = k + 1
+             IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
           ENDDO
           kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 )
           kstart = kflo(kk)
+!AH
+          write(9,"('init kflo, kfuo: ', 4(i3,2x), e12.5)") kk, kflo(kk), kfuo(kk), nzt,  cg%zw(kk)
+          flush(9)
        ENDDO
 !
-!--    Precomputation of number of fine-grid nodes inside coarse-grid ij-faces.
-!--    Note that ii, jj, and kk are coarse-grid indices. 
+!--    Precomputation of number of fine-grid nodes inside parent-grid cells.
+!--    Note that ii, jj, and kk are parent-grid indices. 
 !--    This information is needed in anterpolation.
        DO  ii = icl, icr
           DO  jj = jcs, jcn
-             DO kk = 0, kctu
+!AH             DO kk = 0, kcto+1
+             DO kk = 0, cg%nz+1
 !
 !--             u-component
@@ -3223,5 +3411,6 @@
              ENDDO
 
-             DO kk = 0, kctw
+!AH             DO kk = 0, kctw+1
+             DO kk = 0, cg%nz+1
 !
 !--             w-component 
@@ -3259,5 +3448,4 @@
              ENDIF
           ENDDO
-
 
           DO  jj = jcs, jcn
@@ -3275,6 +3463,6 @@
        ENDIF
       
-       ALLOCATE( fraz(0:kctu) )
-       DO  kk = 0, kctu
+       ALLOCATE( fraz(0:kcto) )
+       DO  kk = 0, kcto
           zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4
           fraz(kk) = zeta / ( 1.0_wp + zeta )
@@ -3462,5 +3650,5 @@
        coord_x(i) = lower_left_coord_x + i * dx
     ENDDO
-     
+
     DO  j = -nbgp, ny + nbgp
        coord_y(j) = lower_left_coord_y + j * dy
@@ -3472,5 +3660,4 @@
 
 
-
  SUBROUTINE pmci_set_array_pointer( name, child_id, nz_cl, n )
 
@@ -3487,4 +3674,5 @@
 
     REAL(wp), POINTER, DIMENSION(:,:)     ::  p_2d        !<
+    REAL(wp), POINTER, DIMENSION(:,:)     ::  p_2d_sec    !<
     REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d        !<
     REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d_sec    !<
@@ -3584,4 +3772,5 @@
 
 
+
 INTEGER FUNCTION get_number_of_childs ()
 
@@ -3614,4 +3803,5 @@
 
 END FUNCTION get_childid
+
 
 
@@ -3641,4 +3831,6 @@
 END SUBROUTINE get_child_edges
 
+
+
 SUBROUTINE  get_child_gridspacing (m, dx,dy,dz)
 
@@ -3656,5 +3848,7 @@
 END SUBROUTINE get_child_gridspacing
 
-SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc,n  )
+
+
+SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc, n  )
 
     IMPLICIT NONE
@@ -3666,5 +3860,5 @@
     INTEGER(iwp), INTENT(IN) ::  je      !<
     INTEGER(iwp), INTENT(IN) ::  js      !<
-    INTEGER(iwp), INTENT(IN) ::  nzc     !<  Note that nzc is cg%nz
+    INTEGER(iwp), INTENT(IN) ::  nzc     !<  nzc is cg%nz, but note that cg%nz is not the original nz of parent, but the highest parent-grid level needed for nesting.
 
     INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species
@@ -3819,9 +4013,12 @@
 !--    The interpolation.
        CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
-                                   r2yo, r1zo, r2zo, 'u' )
+                                   r2yo, r1zo, r2zo, kcto, iflu, ifuu,         &
+                                   jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' )
        CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,   &
-                                   r2yv, r1zo, r2zo, 'v' )
+                                   r2yv, r1zo, r2zo, kcto, iflo, ifuo,         &
+                                   jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' )
        CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,   &
-                                   r2yo, r1zw, r2zw, 'w' )
+                                   r2yo, r1zw, r2zw, kctw, iflo, ifuo,         &
+                                   jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' )
 
        IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.          &
@@ -3829,15 +4026,18 @@
                .NOT. constant_diffusion ) )  THEN
           CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
-                                      r2yo, r1zo, r2zo, 'e' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,       &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' )
        ENDIF
 
        IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
           CALL pmci_interp_tril_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
-                                      r1yo, r2yo, r1zo, r2zo, 's' )
+                                      r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,&
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
        IF ( .NOT. neutral )  THEN
           CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
-                                      r1yo, r2yo, r1zo, r2zo, 's' )
+                                      r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,&
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
@@ -3845,18 +4045,27 @@
 
           CALL pmci_interp_tril_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, &
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
           IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
              CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,   &
-                                          r1yo, r2yo, r1zo, r2zo, 's' ) 
+                                         r1yo, r2yo, r1zo, r2zo, kcto,         &
+                                         iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
+                                         ijkfc_s, 's' ) 
              CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,   &
-                                         r1yo, r2yo, r1zo, r2zo, 's' )   
+                                         r1yo, r2yo, r1zo, r2zo, kcto,         &
+                                         iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
+                                         ijkfc_s, 's' )   
           ENDIF
 
           IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
              CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,   &
-                                         r1yo, r2yo, r1zo, r2zo, 's' )
+                                         r1yo, r2yo, r1zo, r2zo, kcto,         &
+                                         iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
+                                         ijkfc_s, 's' )
              CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,   &
-                                         r1yo, r2yo, r1zo, r2zo, 's' )
+                                         r1yo, r2yo, r1zo, r2zo, kcto,         &
+                                         iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
+                                         ijkfc_s, 's' )
           ENDIF
 
@@ -3865,5 +4074,6 @@
        IF ( passive_scalar )  THEN
           CALL pmci_interp_tril_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,  &
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
@@ -3873,5 +4083,6 @@
                                          chem_spec_c(:,:,:,n),                 &
                                          ico, jco, kco, r1xo, r2xo, r1yo,      &
-                                         r2yo, r1zo, r2zo, 's' )
+                                         r2yo, r1zo, r2zo, kcto, iflo, ifuo,   &
+                                         jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
           ENDDO
        ENDIF
@@ -3912,5 +4123,6 @@
 
     SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
-                                     r1z, r2z, var )
+                                     r1z, r2z, kct, ifl, ifu, jfl, jfu,        &
+                                     kfl, kfu, ijkfc, var )
 !
 !--    Interpolation of the internal values for the child-domain initialization
@@ -3925,7 +4137,26 @@
        INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
 
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !<
+       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc       !<
+       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x   !<
+       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x   !<
+       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y   !<
+       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y   !<
+       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z   !<
+       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z   !<
+
+       INTEGER(iwp) :: kct
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+
        INTEGER(iwp) ::  i        !<
        INTEGER(iwp) ::  ib       !<
        INTEGER(iwp) ::  ie       !<
+       INTEGER(iwp) ::  ijk      !<
        INTEGER(iwp) ::  j        !<
        INTEGER(iwp) ::  jb       !<
@@ -3939,14 +4170,8 @@
        INTEGER(iwp) ::  m        !<
        INTEGER(iwp) ::  n        !<
-
-       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !<
-       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc       !<
-       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x   !<
-       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x   !<
-       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y   !<
-       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y   !<
-       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z   !<
-       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z   !<
-
+       INTEGER(iwp) ::  var_flag !<
+
+       REAL(wp) ::  cellsum    !<
+       REAL(wp) ::  cellsumd   !<
        REAL(wp) ::  fk         !<
        REAL(wp) ::  fkj        !<
@@ -3957,4 +4182,6 @@
        REAL(wp) ::  logratio   !<
        REAL(wp) ::  logzuc1    !<
+       REAL(wp) ::  rcorr      !<
+       REAL(wp) ::  rcorr_ijk  !<       
        REAL(wp) ::  zuc1       !<
        REAL(wp) ::  z0_topo    !<  roughness at vertical walls
@@ -3965,4 +4192,5 @@
        jb = nys
        je = nyn
+       kb = 0
        IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
@@ -3989,4 +4217,14 @@
           ENDIF
        ENDIF
+
+       IF ( var == 'u' )  THEN 
+          var_flag = 1 
+       ELSEIF ( var == 'v' )  THEN
+          var_flag = 2 
+       ELSEIF ( var == 'w' )  THEN
+          var_flag = 3
+       ELSE
+          var_flag = 0
+       ENDIF
 !
 !--    Trilinear interpolation.
@@ -3997,5 +4235,5 @@
 !--          topography top at grid point (j,i) in order to not overwrite 
 !--          the bottom BC.
-             kb = get_topography_top_index_ji( j, i, TRIM ( var ) ) + 1
+!             kb = get_topography_top_index_ji( j, i, TRIM ( var ) ) + 1             
              DO  k = kb, nzt + 1
                 l = ic(i)
@@ -4062,4 +4300,64 @@
 
        ENDIF
+!
+!--    Apply the reversibility correction. 
+       DO  l = icl, icr
+          DO  m = jcs, jcn
+             DO  n = 0, kct+1
+                ijk = 1
+                cellsum   = 0.0_wp
+                cellsumd  = 0.0_wp
+!
+!--             Note that the index name i must not be used here as a loop 
+!--             index name since i is the constant boundary index, hence 
+!--             the name ia. 
+                DO  i = ifl(l), ifu(l)   
+                   DO  j = jfl(m), jfu(m)
+                      DO  k = kfl(n), kfu(n)
+                         cellsum = cellsum + MERGE( f(k,j,i), 0.0_wp,           &
+                              BTEST( wall_flags_0(k,j,i), var_flag ) )
+                         celltmpd(ijk) = ABS( fc(n,m,l) - f(k,j,i) )
+                         cellsumd      = cellsumd  + MERGE( celltmpd(ijk),      &
+                              0.0_wp, BTEST( wall_flags_0(k,j,i), var_flag ) )
+                         ijk = ijk + 1
+                      ENDDO
+                   ENDDO
+                ENDDO
+                
+                IF ( ijkfc(n,m,l) /= 0 )  THEN
+                   cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+                   rcorr     = fc(n,m,l) - cellsum
+                   cellsumd  = cellsumd / REAL( ijkfc(n,m,l), KIND=wp )
+                ELSE
+                   cellsum   = 0.0_wp                  
+                   rcorr     = 0.0_wp
+                   cellsumd  = 1.0_wp
+                   celltmpd  = 1.0_wp
+                ENDIF
+! 
+!--             Distribute the correction term to the child nodes according to 
+!--             their relative difference to the parent value such that the
+!--             node with the largest difference gets the largest share of the 
+!--             correction. The distribution is skipped if rcorr is negligibly 
+!--             small in order to avoid division by zero. 
+                IF ( ABS(rcorr) < 0.000001_wp )  THEN                  
+                   cellsumd  = 1.0_wp
+                   celltmpd  = 1.0_wp
+                ENDIF
+
+                ijk = 1
+                DO  i = ifl(l), ifu(l)
+                   DO  j = jfl(m), jfu(m)
+                      DO  k = kfl(n), kfu(n)
+                         rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+                         f(k,j,i) = f(k,j,i) + rcorr_ijk
+                         ijk = ijk + 1
+                      ENDDO
+                   ENDDO
+                ENDDO
+             
+             ENDDO  ! n
+          ENDDO  ! m
+       ENDDO  ! l
 
     END SUBROUTINE pmci_interp_tril_all
@@ -4165,130 +4463,4 @@
 
 
- SUBROUTINE pmci_ensure_nest_mass_conservation
-
-!
-!-- Adjust the volume-flow rate through the top boundary so that the net volume
-!-- flow through all boundaries of the current nest domain becomes zero.
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  i                           !<
-    INTEGER(iwp) ::  ierr                        !<
-    INTEGER(iwp) ::  j                           !<
-    INTEGER(iwp) ::  k                           !<
-
-    REAL(wp) ::  dxdy                            !<
-    REAL(wp) ::  innor                           !<
-    REAL(wp) ::  w_lt                            !<
-    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !<
-
-!
-!-- Sum up the volume flow through the left/right boundaries
-    volume_flow(1)   = 0.0_wp
-    volume_flow_l(1) = 0.0_wp
-
-    IF ( bc_dirichlet_l )  THEN
-       i = 0
-       innor = dy
-       DO   j = nys, nyn
-          DO   k = nzb+1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-          ENDDO
-       ENDDO
-    ENDIF
-
-    IF ( bc_dirichlet_r )  THEN
-       i = nx + 1
-       innor = -dy
-       DO   j = nys, nyn
-          DO   k = nzb+1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 1 ) )
-          ENDDO
-       ENDDO
-    ENDIF
-
-#if defined( __parallel )
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL,         &
-                        MPI_SUM, comm2d, ierr )
-#else
-    volume_flow(1) = volume_flow_l(1)
-#endif
-    
-!
-!-- Sum up the volume flow through the south/north boundaries
-    volume_flow(2)   = 0.0_wp
-    volume_flow_l(2) = 0.0_wp
-
-    IF ( bc_dirichlet_s )  THEN
-       j = 0
-       innor = dx
-       DO   i = nxl, nxr
-          DO   k = nzb+1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 2 ) )
-          ENDDO
-       ENDDO
-    ENDIF
-
-    IF ( bc_dirichlet_n )  THEN
-       j = ny + 1
-       innor = -dx
-       DO   i = nxl, nxr
-          DO   k = nzb+1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
-                                     * MERGE( 1.0_wp, 0.0_wp,                  &
-                                              BTEST( wall_flags_0(k,j,i), 2 ) )
-          ENDDO
-       ENDDO
-    ENDIF
-
-#if defined( __parallel )
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL,         &
-                        MPI_SUM, comm2d, ierr )
-#else
-    volume_flow(2) = volume_flow_l(2)
-#endif
-
-!
-!-- Sum up the volume flow through the top boundary
-    volume_flow(3)   = 0.0_wp
-    volume_flow_l(3) = 0.0_wp
-    dxdy = dx * dy
-    k = nzt
-    DO   i = nxl, nxr
-       DO   j = nys, nyn
-          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
-       ENDDO
-    ENDDO
-
-#if defined( __parallel )
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL,         &
-                        MPI_SUM, comm2d, ierr )
-#else
-    volume_flow(3) = volume_flow_l(3)
-#endif
-
-!
-!-- Correct the top-boundary value of w
-    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
-    DO   i = nxl, nxr
-       DO   j = nys, nyn
-          DO  k = nzt, nzt + 1
-             w(k,j,i) = w(k,j,i) + w_lt
-          ENDDO
-       ENDDO
-    ENDDO
-
- END SUBROUTINE pmci_ensure_nest_mass_conservation
-
-
-
  SUBROUTINE pmci_synchronize
 
@@ -4409,5 +4581,4 @@
     DO  m = 1, SIZE( pmc_parent_for_child ) - 1
        child_id = pmc_parent_for_child(m)
-       
        IF ( direction == parent_to_child )  THEN
           CALL cpu_log( log_point_s(71), 'pmc parent send', 'start' )
@@ -4484,5 +4655,5 @@
 
        IF ( direction == parent_to_child )  THEN
-
+    
           CALL cpu_log( log_point_s(73), 'pmc child recv', 'start' )
           CALL pmc_c_getbuffer( )
@@ -4492,5 +4663,5 @@
           CALL pmci_interpolation
           CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' )
-
+     
        ELSE
 !
@@ -4517,4 +4688,6 @@
        IMPLICIT NONE
 
+       INTEGER(iwp) ::  ibgp       !< index running over the nbgp boundary ghost points in i-direction
+       INTEGER(iwp) ::  jbgp       !< index running over the nbgp boundary ghost points in j-direction
        INTEGER(iwp) ::  n          !< running index for number of chemical species
       
@@ -4527,31 +4700,41 @@
 !--       Left border pe:
           IF ( bc_dirichlet_l )  THEN
-             
-             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
-                                       r1yo, r2yo, r1zo, r2zo,                 &
-                                       logc_u_l, logc_ratio_u_l,               &
-                                       logc_kbounds_u_l,                       &
-                                       nzt_topo_nestbc_l, 'l', 'u' )
-
-             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
-                                       r1yv, r2yv, r1zo, r2zo,                 &
-                                       logc_v_l, logc_ratio_v_l,               &
-                                       logc_kbounds_v_l,                       &               
-                                       nzt_topo_nestbc_l, 'l', 'v' )
-
-             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
-                                       r1yo, r2yo, r1zw, r2zw,                 &
-                                       logc_w_l, logc_ratio_w_l,               &
-                                       logc_kbounds_w_l,                       &
-                                       nzt_topo_nestbc_l, 'l', 'w' )
+
+             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
+                                       r1yo, r2yo, r1zo, r2zo,                  &
+                                       logc_u_l, logc_ratio_u_l,                &
+                                       logc_kbounds_u_l, nzt_topo_nestbc_l,     &
+                                       kcto, iflu, ifuu, jflo, jfuo, kflo,      &
+                                       kfuo, ijkfc_u, 'l', 'u' )
+
+             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
+                                       r1yv, r2yv, r1zo, r2zo,                  &
+                                       logc_v_l, logc_ratio_v_l,                &
+                                       logc_kbounds_v_l, nzt_topo_nestbc_l,     &
+                                       kcto, iflo, ifuo, jflv, jfuv, kflo,      &
+                                       kfuo, ijkfc_v, 'l', 'v' )
+
+             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
+                                       r1yo, r2yo, r1zw, r2zw,                  &
+                                       logc_w_l, logc_ratio_w_l,                &
+                                       logc_kbounds_w_l, nzt_topo_nestbc_l,     &
+                                       kctw, iflo, ifuo, jflo, jfuo, kflw,      &
+                                       kfuw, ijkfc_w, 'l', 'w' )
 
              IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
                   (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
                      .NOT. constant_diffusion ) )  THEN
-                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo,              &
-                                          logc_w_l, logc_ratio_w_l,            &
-                                          logc_kbounds_w_l,                    &
-                                          nzt_topo_nestbc_l, 'l', 'e' )
+!                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
+!                                          r1yo, r2yo, r1zo, r2zo,              &
+!                                          logc_w_l, logc_ratio_w_l,            &
+!                                          logc_kbounds_w_l, nzt_topo_nestbc_l, &
+!                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+!                                          kfuo, ijkfc_s, 'l', 'e' )
+!
+!--             Interpolation of e is replaced by the Neumann condition. 
+                DO ibgp = -nbgp, -1
+                   e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,0)
+                ENDDO
+
              ENDIF
 
@@ -4560,6 +4743,7 @@
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_l, logc_ratio_w_l,            &
-                                          logc_kbounds_w_l,                    &
-                                          nzt_topo_nestbc_l, 'l', 's' )
+                                          logc_kbounds_w_l, nzt_topo_nestbc_l, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'l', 's' )
              ENDIF
 
@@ -4568,6 +4752,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_l, logc_ratio_w_l,            &
-                                          logc_kbounds_w_l,                    &                
-                                          nzt_topo_nestbc_l, 'l', 's' )
+                                          logc_kbounds_w_l, nzt_topo_nestbc_l, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'l', 's' )
              ENDIF
 
@@ -4577,6 +4762,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_l, logc_ratio_w_l,            &
-                                          logc_kbounds_w_l,                    &
-                                          nzt_topo_nestbc_l, 'l', 's' )
+                                          logc_kbounds_w_l, nzt_topo_nestbc_l, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'l', 's' )
 
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
@@ -4585,5 +4771,7 @@
                                              logc_w_l, logc_ratio_w_l,         &
                                              logc_kbounds_w_l,                 &
-                                             nzt_topo_nestbc_l, 'l', 's' )  
+                                             nzt_topo_nestbc_l,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'l', 's' )  
 
                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
@@ -4591,5 +4779,7 @@
                                              logc_w_l, logc_ratio_w_l,         &
                                              logc_kbounds_w_l,                 &
-                                             nzt_topo_nestbc_l, 'l', 's' )          
+                                             nzt_topo_nestbc_l,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'l', 's' )          
                 ENDIF
 
@@ -4599,11 +4789,15 @@
                                              logc_w_l, logc_ratio_w_l,         &
                                              logc_kbounds_w_l,                 &
-                                             nzt_topo_nestbc_l, 'l', 's' ) 
+                                             nzt_topo_nestbc_l,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'l', 's' ) 
 
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_l, logc_ratio_w_l,         &
-                                             logc_kbounds_w_l,                 &                
-                                             nzt_topo_nestbc_l, 'l', 's' )             
+                                             logc_kbounds_w_l,                 &
+                                             nzt_topo_nestbc_l,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'l', 's' )             
                 ENDIF
 
@@ -4615,5 +4809,7 @@
                                           logc_w_l, logc_ratio_w_l,            &
                                           logc_kbounds_w_l,                    & 
-                                          nzt_topo_nestbc_l, 'l', 's' )
+                                          nzt_topo_nestbc_l,                   &
+                                          kcto, iflo, ifuo, jflo, jfuo,        &
+                                          kflo, kfuo, ijkfc_s, 'l', 's' )
              ENDIF
 
@@ -4626,5 +4822,7 @@
                                              logc_w_l, logc_ratio_w_l,         &
                                              logc_kbounds_w_l,                 & 
-                                             nzt_topo_nestbc_l, 'l', 's' )
+                                             nzt_topo_nestbc_l,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'l', 's' )
                 ENDDO
              ENDIF
@@ -4638,28 +4836,36 @@
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_r, logc_ratio_u_r,               &
-                                       logc_kbounds_u_r,                       &
-                                       nzt_topo_nestbc_r, 'r', 'u' )
+                                       logc_kbounds_u_r, nzt_topo_nestbc_r,    &
+                                       kcto, iflu, ifuu, jflo, jfuo, kflo,     &
+                                       kfuo, ijkfc_u, 'r', 'u' )
 
              CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
                                        r1yv, r2yv, r1zo, r2zo,                 &
                                        logc_v_r, logc_ratio_v_r,               &
-                                       logc_kbounds_v_r,                       &
-                                       nzt_topo_nestbc_r, 'r', 'v' )
+                                       logc_kbounds_v_r, nzt_topo_nestbc_r,    &
+                                       kcto, iflo, ifuo, jflv, jfuv, kflo,     &
+                                       kfuo, ijkfc_v, 'r', 'v' )
 
              CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
                                        r1yo, r2yo, r1zw, r2zw,                 &
                                        logc_w_r, logc_ratio_w_r,               &
-                                       logc_kbounds_w_r,                       &
-                                       nzt_topo_nestbc_r, 'r', 'w' )
+                                       logc_kbounds_w_r, nzt_topo_nestbc_r,    &
+                                       kctw, iflo, ifuo, jflo, jfuo, kflw,     &
+                                       kfuw, ijkfc_w, 'r', 'w' )
 
              IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
                   (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
                      .NOT. constant_diffusion ) )  THEN
-                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo,r2yo, r1zo, r2zo,               &
-                                          logc_w_r, logc_ratio_w_r,            &
-                                          logc_kbounds_w_r,                    &
-                                          nzt_topo_nestbc_r, 'r', 'e' )
-
+!                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
+!                                          r1yo,r2yo, r1zo, r2zo,               &
+!                                          logc_w_r, logc_ratio_w_r,            &
+!                                          logc_kbounds_w_r, nzt_topo_nestbc_r, &
+!                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+!                                          kfuo, ijkfc_s, 'r', 'e' )
+!
+!--             Interpolation of e is replaced by the Neumann condition. 
+                DO ibgp = nx+1, nx+nbgp
+                   e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,nx)
+                ENDDO
              ENDIF
 
@@ -4668,6 +4874,7 @@
                                           r2xo, r1yo,r2yo, r1zo, r2zo,         &
                                           logc_w_r, logc_ratio_w_r,            &
-                                          logc_kbounds_w_r,                    &
-                                          nzt_topo_nestbc_r, 'r', 's' )
+                                          logc_kbounds_w_r, nzt_topo_nestbc_r, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'r', 's' )
 
              ENDIF
@@ -4677,6 +4884,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_r, logc_ratio_w_r,            &
-                                          logc_kbounds_w_r,                    &
-                                          nzt_topo_nestbc_r, 'r', 's' )
+                                          logc_kbounds_w_r, nzt_topo_nestbc_r, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'r', 's' )
              ENDIF
 
@@ -4685,6 +4893,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_r, logc_ratio_w_r,            &
-                                          logc_kbounds_w_r,                    &
-                                          nzt_topo_nestbc_r, 'r', 's' )
+                                          logc_kbounds_w_r, nzt_topo_nestbc_r, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'r', 's' )
 
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
@@ -4694,5 +4903,7 @@
                                              logc_w_r, logc_ratio_w_r,         &
                                              logc_kbounds_w_r,                 &
-                                             nzt_topo_nestbc_r, 'r', 's' ) 
+                                             nzt_topo_nestbc_r,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'r', 's' ) 
      
                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
@@ -4700,6 +4911,7 @@
                                              logc_w_r, logc_ratio_w_r,         &
                                              logc_kbounds_w_r,                 &
-                                             nzt_topo_nestbc_r, 'r', 's' )
-
+                                             nzt_topo_nestbc_r,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'r', 's' )
 
                 ENDIF
@@ -4712,5 +4924,8 @@
                                              logc_w_r, logc_ratio_w_r,         &
                                              logc_kbounds_w_r,                 &
-                                             nzt_topo_nestbc_r, 'r', 's' ) 
+                                             nzt_topo_nestbc_r,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s,              &
+                                             'r', 's' ) 
 
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
@@ -4718,5 +4933,7 @@
                                              logc_w_r, logc_ratio_w_r,         &
                                              logc_kbounds_w_r,                 &
-                                             nzt_topo_nestbc_r, 'r', 's' ) 
+                                             nzt_topo_nestbc_r,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'r', 's' ) 
 
                 ENDIF
@@ -4729,5 +4946,7 @@
                                           logc_w_r, logc_ratio_w_r,            &
                                           logc_kbounds_w_r,                    &
-                                          nzt_topo_nestbc_r, 'r', 's' )
+                                          nzt_topo_nestbc_r,                   &
+                                          kcto, iflo, ifuo, jflo, jfuo,        &
+                                          kflo, kfuo, ijkfc_s, 'r', 's' )
              ENDIF
 
@@ -4740,5 +4959,7 @@
                                              logc_w_r, logc_ratio_w_r,         &
                                              logc_kbounds_w_r,                 &
-                                             nzt_topo_nestbc_r, 'r', 's' )
+                                             nzt_topo_nestbc_r,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'r', 's' )
                 ENDDO
              ENDIF
@@ -4751,28 +4972,36 @@
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_s, logc_ratio_u_s,               &
-                                       logc_kbounds_u_s,                       &
-                                       nzt_topo_nestbc_s, 's', 'u' )
+                                       logc_kbounds_u_s, nzt_topo_nestbc_s,    &
+                                       kcto, iflu, ifuu, jflo, jfuo, kflo,     &
+                                       kfuo, ijkfc_u, 's', 'u' )
 
              CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
                                        r1yv, r2yv, r1zo, r2zo,                 &
                                        logc_v_s, logc_ratio_v_s,               &
-                                       logc_kbounds_v_s,                       &
-                                       nzt_topo_nestbc_s, 's', 'v' )
+                                       logc_kbounds_v_s, nzt_topo_nestbc_s,    &
+                                       kcto, iflo, ifuo, jflv, jfuv, kflo,     &
+                                       kfuo, ijkfc_v, 's', 'v' )
 
              CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
                                        r1yo, r2yo, r1zw, r2zw,                 &
                                        logc_w_s, logc_ratio_w_s,               &
-                                       logc_kbounds_w_s,                       &
-                                       nzt_topo_nestbc_s, 's','w' )
+                                       logc_kbounds_w_s, nzt_topo_nestbc_s,    &
+                                       kctw, iflo, ifuo, jflo, jfuo, kflw,     &
+                                       kfuw, ijkfc_w, 's','w' )
 
              IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
                   (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
                      .NOT. constant_diffusion ) )  THEN
-                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo,              &
-                                          logc_w_s, logc_ratio_w_s,            &
-                                          logc_kbounds_w_s,                    &
-                                          nzt_topo_nestbc_s, 's', 'e' )
-
+!                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
+!                                          r1yo, r2yo, r1zo, r2zo,              &
+!                                          logc_w_s, logc_ratio_w_s,            &
+!                                          logc_kbounds_w_s, nzt_topo_nestbc_s, &
+!                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+!                                          kfuo, ijkfc_s, 's', 'e' )
+!
+!--             Interpolation of e is replaced by the Neumann condition. 
+                DO jbgp = -nbgp, -1
+                   e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,0,nxl:nxr)
+                ENDDO
              ENDIF
 
@@ -4781,7 +5010,7 @@
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_s, logc_ratio_w_s,            &
-                                          logc_kbounds_w_s,                    &
-                                          nzt_topo_nestbc_s, 's', 's' )
-
+                                          logc_kbounds_w_s, nzt_topo_nestbc_s, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 's', 's' )
              ENDIF
 
@@ -4790,6 +5019,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_s, logc_ratio_w_s,            &
-                                          logc_kbounds_w_s,                    &
-                                          nzt_topo_nestbc_s, 's', 's' )
+                                          logc_kbounds_w_s, nzt_topo_nestbc_s, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 's', 's' )
              ENDIF
 
@@ -4798,6 +5028,7 @@
                                           r1yo,r2yo, r1zo, r2zo,               &
                                           logc_w_s, logc_ratio_w_s,            &
-                                          logc_kbounds_w_s,                    &
-                                          nzt_topo_nestbc_s, 's', 's' )
+                                          logc_kbounds_w_s, nzt_topo_nestbc_s, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 's', 's' )
 
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
@@ -4807,5 +5038,7 @@
                                              logc_w_s, logc_ratio_w_s,         &
                                              logc_kbounds_w_s,                 &
-                                             nzt_topo_nestbc_s, 's', 's' )
+                                             nzt_topo_nestbc_s,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 's', 's' )
 
                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
@@ -4813,5 +5046,7 @@
                                              logc_w_s, logc_ratio_w_s,         &
                                              logc_kbounds_w_s,                 &
-                                             nzt_topo_nestbc_s, 's', 's' )
+                                             nzt_topo_nestbc_s,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 's', 's' )
 
                 ENDIF
@@ -4823,5 +5058,7 @@
                                              logc_w_s, logc_ratio_w_s,         &
                                              logc_kbounds_w_s,                 &
-                                             nzt_topo_nestbc_s, 's', 's' )
+                                             nzt_topo_nestbc_s,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 's', 's' )
 
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
@@ -4829,5 +5066,7 @@
                                              logc_w_s, logc_ratio_w_s,         &
                                              logc_kbounds_w_s,                 &
-                                             nzt_topo_nestbc_s, 's', 's' )
+                                             nzt_topo_nestbc_s,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 's', 's' )
 
                 ENDIF
@@ -4840,5 +5079,7 @@
                                           logc_w_s, logc_ratio_w_s,            &
                                           logc_kbounds_w_s,                    &
-                                          nzt_topo_nestbc_s, 's', 's' )
+                                          nzt_topo_nestbc_s,                   &
+                                          kcto, iflo, ifuo, jflo, jfuo,        &
+                                          kflo, kfuo, ijkfc_s, 's', 's' )
              ENDIF
 
@@ -4851,5 +5092,7 @@
                                              logc_w_s, logc_ratio_w_s,         &
                                              logc_kbounds_w_s,                 &
-                                             nzt_topo_nestbc_s, 's', 's' )
+                                             nzt_topo_nestbc_s,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 's', 's' )
                 ENDDO
              ENDIF
@@ -4862,28 +5105,36 @@
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_n, logc_ratio_u_n,               &
-                                       logc_kbounds_u_n,                       &
-                                       nzt_topo_nestbc_n, 'n', 'u' )
+                                       logc_kbounds_u_n, nzt_topo_nestbc_n,    &
+                                       kcto, iflu, ifuu, jflo, jfuo, kflo,     &
+                                       kfuo, ijkfc_u, 'n', 'u' )
 
              CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
                                        r1yv, r2yv, r1zo, r2zo,                 &
                                        logc_v_n, logc_ratio_v_n,               &
-                                       logc_kbounds_v_n,                       &
-                                       nzt_topo_nestbc_n, 'n', 'v' )
+                                       logc_kbounds_v_n, nzt_topo_nestbc_n,    &
+                                       kcto, iflo, ifuo, jflv, jfuv, kflo,     &
+                                       kfuo, ijkfc_v, 'n', 'v' )
 
              CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
                                        r1yo, r2yo, r1zw, r2zw,                 &
                                        logc_w_n, logc_ratio_w_n,               &
-                                       logc_kbounds_w_n,                       &
-                                       nzt_topo_nestbc_n, 'n', 'w' )
+                                       logc_kbounds_w_n, nzt_topo_nestbc_n,    &
+                                       kctw, iflo, ifuo, jflo, jfuo, kflw,     &
+                                       kfuw, ijkfc_w, 'n', 'w' )
 
              IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   & 
                   (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
                      .NOT. constant_diffusion ) )  THEN
-                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
-                                          r1yo, r2yo, r1zo, r2zo,              &
-                                          logc_w_n, logc_ratio_w_n,            &
-                                          logc_kbounds_w_n,                    &
-                                          nzt_topo_nestbc_n, 'n', 'e' )
-
+!                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
+!                                          r1yo, r2yo, r1zo, r2zo,              &
+!                                          logc_w_n, logc_ratio_w_n,            &
+!                                          logc_kbounds_w_n, nzt_topo_nestbc_n, &
+!                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+!                                          kfuo, ijkfc_s, 'n', 'e' )
+!
+!--             Interpolation of e is replaced by the Neumann condition. 
+                DO jbgp = ny+1, ny+nbgp
+                   e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,ny,nxl:nxr)
+                ENDDO
              ENDIF
 
@@ -4892,6 +5143,7 @@
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_n, logc_ratio_w_n,            &
-                                          logc_kbounds_w_n,                    &
-                                          nzt_topo_nestbc_n, 'n', 's' )
+                                          logc_kbounds_w_n, nzt_topo_nestbc_n, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'n', 's' )
 
              ENDIF
@@ -4901,6 +5153,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_n, logc_ratio_w_n,            &
-                                          logc_kbounds_w_n,                    &
-                                          nzt_topo_nestbc_n, 'n', 's' )
+                                          logc_kbounds_w_n, nzt_topo_nestbc_n, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'n', 's' )
              ENDIF
 
@@ -4909,6 +5162,7 @@
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_n, logc_ratio_w_n,            &
-                                          logc_kbounds_w_n,                    &
-                                          nzt_topo_nestbc_n, 'n', 's' )
+                                          logc_kbounds_w_n, nzt_topo_nestbc_n, &
+                                          kcto, iflo, ifuo, jflo, jfuo, kflo,  &
+                                          kfuo, ijkfc_s, 'n', 's' )
 
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
@@ -4918,5 +5172,7 @@
                                              logc_w_n, logc_ratio_w_n,         &
                                              logc_kbounds_w_n,                 &
-                                             nzt_topo_nestbc_n, 'n', 's' )
+                                             nzt_topo_nestbc_n,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'n', 's' )
 
                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
@@ -4924,5 +5180,7 @@
                                              logc_u_n, logc_ratio_u_n,         &
                                              logc_kbounds_w_n,                 &
-                                             nzt_topo_nestbc_n, 'n', 's' )
+                                             nzt_topo_nestbc_n,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'n', 's' )
 
                 ENDIF
@@ -4934,5 +5192,7 @@
                                              logc_w_n, logc_ratio_w_n,         &
                                              logc_kbounds_w_n,                 &
-                                             nzt_topo_nestbc_n, 'n', 's' )
+                                             nzt_topo_nestbc_n,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'n', 's' )
 
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
@@ -4940,5 +5200,7 @@
                                              logc_w_n, logc_ratio_w_n,         &
                                              logc_kbounds_w_n,                 &
-                                             nzt_topo_nestbc_n, 'n', 's' )
+                                             nzt_topo_nestbc_n,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'n', 's' )
 
                 ENDIF
@@ -4951,5 +5213,7 @@
                                           logc_w_n, logc_ratio_w_n,            &
                                           logc_kbounds_w_n,                    &
-                                          nzt_topo_nestbc_n, 'n', 's' )
+                                          nzt_topo_nestbc_n,                   &
+                                          kcto, iflo, ifuo, jflo, jfuo,        &
+                                          kflo, kfuo, ijkfc_s, 'n', 's' )
              ENDIF
 
@@ -4962,5 +5226,7 @@
                                              logc_w_n, logc_ratio_w_n,         &
                                              logc_kbounds_w_n,                 &
-                                             nzt_topo_nestbc_n, 'n', 's' )
+                                             nzt_topo_nestbc_n,                &
+                                             kcto, iflo, ifuo, jflo, jfuo,     &
+                                             kflo, kfuo, ijkfc_s, 'n', 's' )
                 ENDDO
              ENDIF
@@ -4971,25 +5237,35 @@
 !--    All PEs are top-border PEs
        CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
-                                r2yo, r1zo, r2zo, 'u' )
+                                r2yo, r1zo, r2zo, kcto, iflu, ifuu,            &
+                                jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' )
        CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,      &
-                                r2yv, r1zo, r2zo, 'v' )
+                                r2yv, r1zo, r2zo, kcto, iflo, ifuo,            &
+                                jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' )
        CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
-                                r2yo, r1zw, r2zw, 'w' )
+                                r2yo, r1zw, r2zw, kctw, iflo, ifuo,            &
+                                jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' )
 
        IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
             (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
                .NOT. constant_diffusion ) )  THEN
-          CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
-                                   r2yo, r1zo, r2zo, 'e' )
+!          CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
+!                                   r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
+!                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' )
+!
+!--       Interpolation of e is replaced by the Neumann condition. 
+          e(nzt+1,nys:nyn,nxl:nxr) = e(nzt,nys:nyn,nxl:nxr)
+
        ENDIF
 
        IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
           CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
-                                   r1yo, r2yo, r1zo, r2zo, 's' )
+                                   r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,   &
+                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
        IF ( .NOT. neutral )  THEN
-          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
-                                   r2yo, r1zo, r2zo, 's' )
+          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo,         &
+                                   r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,   &
+                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
@@ -4997,13 +5273,16 @@
 
           CALL pmci_interp_tril_t( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo,    &
-                                   r2yo, r1zo, r2zo, 's' )
+                                   r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
+                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
           IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
              CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
              CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
           ENDIF
@@ -5013,8 +5292,10 @@
 
              CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
              CALL pmci_interp_tril_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,&
-                                      r2yo, r1zo, r2zo, 's' )
+                                      r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
+                                      jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
 
           ENDIF
@@ -5024,5 +5305,6 @@
        IF ( passive_scalar )  THEN
           CALL pmci_interp_tril_t( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,     &
-                                   r2yo, r1zo, r2zo, 's' )
+                                   r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
+                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
        ENDIF
 
@@ -5033,5 +5315,6 @@
                                       ico, jco, kco, r1xo, r2xo,               &
                                       r1yo, r2yo, r1zo, r2zo,                  &
-                                      's' )
+                                      kcto, iflo, ifuo, jflo, jfuo,            &
+                                      kflo, kfuo, ijkfc_s, 's' )
           ENDDO
        ENDIF
@@ -5050,9 +5333,7 @@
       INTEGER(iwp) ::  n          !< running index for number of chemical species
 
-
-
-      CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
+      CALL pmci_anterp_tophat( u,  uc,  kcto, iflu, ifuu, jflo, jfuo, kflo,    &
                                kfuo, ijkfc_u, 'u' )
-      CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
+      CALL pmci_anterp_tophat( v,  vc,  kcto, iflo, ifuo, jflv, jfuv, kflo,    &
                                kfuo, ijkfc_v, 'v' )
       CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
@@ -5062,10 +5343,10 @@
 !--   RANS mode.
       IF ( rans_mode_parent  .AND.  rans_mode )  THEN
-         CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
+         CALL pmci_anterp_tophat( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo,   &
                                   kfuo, ijkfc_s, 'e' )
 !
 !--      Anterpolation of dissipation rate only if TKE-e closure is applied.
          IF ( rans_tke_e )  THEN
-            CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
+            CALL pmci_anterp_tophat( diss, dissc, kcto, iflo, ifuo, jflo, jfuo,&
                                      kflo, kfuo, ijkfc_s, 'diss' )
          ENDIF
@@ -5074,5 +5355,5 @@
 
       IF ( .NOT. neutral )  THEN
-         CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo, &
+         CALL pmci_anterp_tophat( pt, ptc, kcto, iflo, ifuo, jflo, jfuo, kflo, &
                                   kfuo, ijkfc_s, 'pt' )
       ENDIF
@@ -5080,13 +5361,13 @@
       IF ( humidity )  THEN
 
-         CALL pmci_anterp_tophat( q, q_c, kctu, iflo, ifuo, jflo, jfuo, kflo,  &
+         CALL pmci_anterp_tophat( q, q_c, kcto, iflo, ifuo, jflo, jfuo, kflo,  &
                                   kfuo, ijkfc_s, 'q' )
 
          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
-            CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
+            CALL pmci_anterp_tophat( qc, qcc, kcto, iflo, ifuo, jflo, jfuo,    &
                                      kflo, kfuo, ijkfc_s, 'qc' )
 
-            CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
+            CALL pmci_anterp_tophat( nc, ncc, kcto, iflo, ifuo, jflo, jfuo,    &
                                      kflo, kfuo, ijkfc_s, 'nc' )
 
@@ -5095,8 +5376,8 @@
          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
 
-            CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
+            CALL pmci_anterp_tophat( qr, qrc, kcto, iflo, ifuo, jflo, jfuo,    &
                                      kflo, kfuo, ijkfc_s, 'qr' )
 
-            CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &
+            CALL pmci_anterp_tophat( nr, nrc, kcto, iflo, ifuo, jflo, jfuo,    &
                                      kflo, kfuo, ijkfc_s, 'nr' )
 
@@ -5106,5 +5387,5 @@
 
       IF ( passive_scalar )  THEN
-         CALL pmci_anterp_tophat( s, sc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
+         CALL pmci_anterp_tophat( s, sc, kcto, iflo, ifuo, jflo, jfuo, kflo,   &
                                   kfuo, ijkfc_s, 's' )
       ENDIF
@@ -5114,5 +5395,5 @@
             CALL pmci_anterp_tophat( chem_species(n)%conc,                     &
                                      chem_spec_c(:,:,:,n),                     &
-                                     kctu, iflo, ifuo, jflo, jfuo, kflo,       &
+                                     kcto, iflo, ifuo, jflo, jfuo, kflo,       &
                                      kfuo, ijkfc_s, 's' )
          ENDDO
@@ -5125,5 +5406,7 @@
    SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
                                    r2z, logc, logc_ratio, logc_kbounds,        &
-                                   nzt_topo_nestbc, edge, var )
+                                   nzt_topo_nestbc,                            &
+                                   kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc,   &
+                                   edge, var )
 !
 !--   Interpolation of ghost-node values used as the child-domain boundary
@@ -5132,6 +5415,4 @@
 
       IMPLICIT NONE
-
-      INTEGER(iwp) ::  nzt_topo_nestbc   !<
 
       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
@@ -5152,6 +5433,21 @@
       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !<
       INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nys:nyn),                &
-                                          INTENT(IN)           ::  logc   !<
-      INTEGER(iwp), DIMENSION(1:2,nys:nyn), INTENT(IN)         ::  logc_kbounds !<
+                                          INTENT(IN)           :: logc   !<
+      INTEGER(iwp), DIMENSION(1:2,nys:nyn), INTENT(IN)         :: logc_kbounds !<
+
+      INTEGER(iwp) :: kct
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+
+!AH      INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+      INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+
+      INTEGER(iwp) ::  nzt_topo_nestbc   !<
 
       CHARACTER(LEN=1), INTENT(IN) ::  edge   !<
@@ -5159,10 +5455,14 @@
 
       INTEGER(iwp) ::  i        !<
+      INTEGER(iwp) ::  ia       !<
       INTEGER(iwp) ::  ib       !<
       INTEGER(iwp) ::  ibgp     !<
+      INTEGER(iwp) ::  ijk      !<
+      INTEGER(iwp) ::  iw       !<
       INTEGER(iwp) ::  j        !<
       INTEGER(iwp) ::  jco      !<
       INTEGER(iwp) ::  jcorr    !<
       INTEGER(iwp) ::  jinc     !<
+      INTEGER(iwp) ::  jw       !<
       INTEGER(iwp) ::  j1       !<
       INTEGER(iwp) ::  k        !<
@@ -5174,5 +5474,9 @@
       INTEGER(iwp) ::  m        !<
       INTEGER(iwp) ::  n        !<
-      
+      INTEGER(iwp) ::  kbc      !<
+      INTEGER(iwp) ::  var_flag !<      
+
+      REAL(wp) ::  cellsum     !<
+      REAL(wP) ::  cellsumd    !<
       REAL(wp) ::  fkj         !<
       REAL(wp) ::  fkjp        !<
@@ -5181,5 +5485,7 @@
       REAL(wp) ::  fk          !<
       REAL(wp) ::  fkp         !<
-      
+      REAL(wp) ::  rcorr       !<
+      REAL(wp) ::  rcorr_ijk   !<
+ 
 ! 
 !--   Check which edge is to be handled
@@ -5198,11 +5504,21 @@
          ib = nxr + 2
       ENDIF
-      
+
+      IF ( var == 'u' )  THEN 
+         var_flag = 1 
+      ELSEIF ( var == 'v' )  THEN
+         var_flag = 2 
+      ELSEIF ( var == 'w' )  THEN
+         var_flag = 3
+      ELSE
+         var_flag = 0
+      ENDIF
+   
       DO  j = nys, nyn+1
 !
 !--      Determine vertical index of topography top at grid point (j,i)
-         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
-
-         DO  k = k_wall, nzt+1
+!AH         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
+
+         DO  k = nzb, nzt+1 !k_wall, nzt+1
             l = ic(i)
             m = jc(j)
@@ -5235,5 +5551,5 @@
                DO  kcorr = 0, ncorr - 1
                   kco = k + kcorr
-                  f(kco,j,i) = logc_ratio(1,kcorr,k,j) * f(k1,j,i)
+!AH                  f(kco,j,i) = logc_ratio(1,kcorr,k,j) * f(k1,j,i)
                ENDDO
             ENDIF
@@ -5261,5 +5577,5 @@
                         jco = j + jinc * jcorr
                         IF ( jco >= nys .AND. jco <= nyn )  THEN
-                           f(k,jco,i) = logc_ratio(2,jcorr,k,j) * f(k,j1,i)
+!AH                           f(k,jco,i) = logc_ratio(2,jcorr,k,j) * f(k,j1,i)
                         ENDIF
                      ENDDO
@@ -5282,8 +5598,8 @@
                         DO  kcorr = 0, ncorr-1
                            kco = k + kcorr
-                           f(kco,jco,i) = 0.5_wp * ( logc_ratio(1,kcorr,k,j) * &
-                                                     f(k1,j,i)                 &
-                                                   + logc_ratio(2,jcorr,k,j) * &
-                                                     f(k,j1,i) )
+!AH                           f(kco,jco,i) = 0.5_wp * ( logc_ratio(1,kcorr,k,j) * &
+!AH                                                     f(k1,j,i)                 &
+!AH                                                   + logc_ratio(2,jcorr,k,j) * &
+!AH                                                     f(k,j1,i) )
                         ENDDO
                      ENDIF
@@ -5295,31 +5611,71 @@
       ENDIF  ! ( topography /= 'flat' )
 !
-!--   Rescale if f is the TKE.
-      IF ( var == 'e')  THEN
-         IF ( edge == 'l' )  THEN
-            DO  j = nys, nyn + 1
-!
-!--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index_ji( j, i, 's' )
-
-               DO  k = k_wall, nzt + 1
-                  f(k,j,i) = tkefactor_l(k,j) * f(k,j,i)
+!--   Apply the reversibility correction to the boundary-normal velocity-
+!--   component u and the scalars. It must not be applied to the boundary-
+!--   tangential velocity components v and w because their 2-D anterpolation
+!--   cells do not cover all the child-grid nodes on the boundary.
+      IF ( .NOT. ( ( var == 'v' ) .OR. ( var == 'w' ) ) )  THEN
+         l = ic(i)
+         DO  m = jcs, jcn
+            DO  n = 0, kct+1
+               ijk = 1
+               cellsum   = 0.0_wp
+               cellsumd  = 0.0_wp
+!
+!--            Note that the index name i must not be used here as a loop 
+!--            index name since i is the constant boundary index, hence 
+!--            the name ia. 
+               DO  ia = ifl(l), ifu(l)   
+                  DO  j = jfl(m), jfu(m)
+                     DO  k = kfl(n), kfu(n)
+                        cellsum = cellsum + MERGE( f(k,j,ia), 0.0_wp,           &
+                             BTEST( wall_flags_0(k,j,ia), var_flag ) )
+                        celltmpd(ijk) = ABS( fc(n,m,l) - f(k,j,ia) )
+                        cellsumd      = cellsumd  + MERGE( celltmpd(ijk),       &
+                             0.0_wp, BTEST( wall_flags_0(k,j,ia), var_flag ) )
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
                ENDDO
-            ENDDO
-         ELSEIF ( edge == 'r' )  THEN           
-            DO  j = nys, nyn+1
-!
-!--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index_ji( j, i, 's' )
-
-               DO  k = k_wall, nzt+1
-                  f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
+
+               IF ( ijkfc(n,m,l) /= 0 )  THEN
+                  cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+                  rcorr     = fc(n,m,l) - cellsum
+                  cellsumd  = cellsumd / REAL( ijkfc(n,m,l), KIND=wp )
+               ELSE
+                  cellsum   = 0.0_wp                  
+                  rcorr     = 0.0_wp
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+! 
+!--            Distribute the correction term to the child nodes according to 
+!--            their relative difference to the parent value such that the
+!--            node with the largest difference gets the largest share of the 
+!--            correction. The distribution is skipped if rcorr is negligibly 
+!--            small in order to avoid division by zero. 
+               IF ( ABS(rcorr) < 0.000001_wp )  THEN                  
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+
+               ijk = 1
+               DO  ia = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     DO  k = kfl(n), kfu(n)
+                        rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+                        f(k,j,ia) = f(k,j,ia) + rcorr_ijk
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
                ENDDO
-            ENDDO
-         ENDIF
-      ENDIF
+               
+            ENDDO  ! n
+         ENDDO  ! m
+         
+      ENDIF  ! var not v or w 
 !
 !--   Store the boundary values also into the other redundant ghost node layers.
-!--   Please note, in case of only one ghost node layer, e.g. for the PW 
+!--   Note that in case of only one ghost node layer, e.g. for the PW 
 !--   scheme, the following loops will not be entered.
       IF ( edge == 'l' )  THEN
@@ -5339,5 +5695,7 @@
    SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
                                    r2z, logc, logc_ratio, logc_kbounds,        &
-                                   nzt_topo_nestbc, edge, var )
+                                   nzt_topo_nestbc,                            &
+                                   kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc,   &
+                                   edge, var )
 
 !
@@ -5347,6 +5705,4 @@
 
       IMPLICIT NONE
-
-      INTEGER(iwp) ::  nzt_topo_nestbc   !<
 
       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
@@ -5370,4 +5726,18 @@
       INTEGER(iwp), DIMENSION(1:2,nxl:nxr), INTENT(IN)         ::  logc_kbounds  !< 
 
+      INTEGER(iwp) :: kct
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH      INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+      INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+
+      INTEGER(iwp) ::  nzt_topo_nestbc   !<
+
       CHARACTER(LEN=1), INTENT(IN) ::  edge   !<
       CHARACTER(LEN=1), INTENT(IN) ::  var    !<
@@ -5377,10 +5747,12 @@
       INTEGER(iwp) ::  icorr   !<
       INTEGER(iwp) ::  ico     !<
+      INTEGER(iwp) ::  ijk     !<
       INTEGER(iwp) ::  i1      !<
       INTEGER(iwp) ::  j       !<
+      INTEGER(iwp) ::  ja      !<
       INTEGER(iwp) ::  jb      !<
       INTEGER(iwp) ::  jbgp    !<
       INTEGER(iwp) ::  k       !<
-      INTEGER(iwp) ::  k_wall   !< vertical index of topography top
+      INTEGER(iwp) ::  k_wall  !< vertical index of topography top
       INTEGER(iwp) ::  kcorr   !<
       INTEGER(iwp) ::  kco     !<
@@ -5389,5 +5761,8 @@
       INTEGER(iwp) ::  m       !<
       INTEGER(iwp) ::  n       !<
-                            
+      INTEGER(iwp) ::  var_flag !<
+
+      REAL(wp) ::  cellsum     !<
+      REAL(wp) ::  cellsumd    !<
       REAL(wp) ::  fk          !<
       REAL(wp) ::  fkj         !<
@@ -5396,4 +5771,6 @@
       REAL(wp) ::  fkpjp       !<
       REAL(wp) ::  fkp         !<
+      REAL(wp) ::  rcorr       !< 
+      REAL(wp) ::  rcorr_ijk   !< 
       
 !
@@ -5414,10 +5791,20 @@
       ENDIF
 
+      IF ( var == 'u' )  THEN 
+         var_flag = 1 
+      ELSEIF ( var == 'v' )  THEN
+         var_flag = 2 
+      ELSEIF ( var == 'w' )  THEN
+         var_flag = 3
+      ELSE
+         var_flag = 0
+      ENDIF
+
       DO  i = nxl, nxr+1
 !
 !--      Determine vertical index of topography top at grid point (j,i)
-         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
-
-         DO  k = k_wall, nzt+1
+!AH         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
+
+         DO  k = nzb, nzt+1 !AH k_wall, nzt+1
             l = ic(i)
             m = jc(j)
@@ -5450,5 +5837,5 @@
                DO  kcorr = 0, ncorr-1
                   kco = k + kcorr
-                  f(kco,j,i) = logc_ratio(1,kcorr,k,i) * f(k1,j,i)
+!AH                  f(kco,j,i) = logc_ratio(1,kcorr,k,i) * f(k1,j,i)
                ENDDO
             ENDIF
@@ -5477,5 +5864,5 @@
                         ico = i + iinc * icorr
                         IF ( ico >= nxl .AND. ico <= nxr )  THEN
-                           f(k,j,ico) = logc_ratio(2,icorr,k,i) * f(k,j,i1)
+!AH                           f(k,j,ico) = logc_ratio(2,icorr,k,i) * f(k,j,i1)
                         ENDIF
                      ENDDO
@@ -5498,8 +5885,8 @@
                         DO  kcorr = 0, ncorr-1
                            kco = k + kcorr
-                           f(kco,j,ico) = 0.5_wp * ( logc_ratio(1,kcorr,k,i) * &
-                                                     f(k1,j,i)                 &
-                                                   + logc_ratio(2,icorr,k,i) * &
-                                                     f(k,j,i1) )
+!AH                           f(kco,j,ico) = 0.5_wp * ( logc_ratio(1,kcorr,k,i) * &
+!AH                                                     f(k1,j,i)                 &
+!AH                                                   + logc_ratio(2,icorr,k,i) * &
+!AH                                                     f(k,j,i1) )
                         ENDDO
                      ENDIF
@@ -5511,29 +5898,67 @@
       ENDIF  ! ( topography /= 'flat' )
 !
-!--   Rescale if f is the TKE.
-      IF ( var == 'e')  THEN
-         IF ( edge == 's' )  THEN
-            DO  i = nxl, nxr + 1
-!
-!--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index_ji( j, i, 's' )
-               DO  k = k_wall, nzt+1
-                  f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
+!--   Apply the reversibility correction to the boundary-normal velocity-
+!--   component v and the scalars. It must not be applied to the boundary-
+!--   tangential velocity components u and w because their 2-D anterpolation 
+!--   cells do not cover all the child-grid nodes on the boundary.
+      IF ( .NOT. ( ( var == 'u' ) .OR. ( var == 'w' ) ) )  THEN
+         m = jc(j)
+         DO  l = icl, icr
+            DO  n = 0, kct+1
+               ijk = 1
+               cellsum   = 0.0_wp 
+               cellsumd  = 0.0_wp
+               DO  i = ifl(l), ifu(l)
+                  DO  ja = jfl(m), jfu(m)
+                     DO  k = kfl(n), kfu(n)
+                        cellsum = cellsum + MERGE( f(k,ja,i), 0.0_wp,           &
+                             BTEST( wall_flags_0(k,ja,i), var_flag ) )                        
+                        celltmpd(ijk) = ABS( fc(n,m,l) - f(k,ja,i) )
+                        cellsumd      = cellsumd  + MERGE( celltmpd(ijk),       &
+                             0.0_wp, BTEST( wall_flags_0(k,ja,i), var_flag ) )
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
                ENDDO
-            ENDDO
-         ELSEIF ( edge == 'n' )  THEN
-            DO  i = nxl, nxr + 1
-!
-!--            Determine vertical index of topography top at grid point (j,i)
-               k_wall = get_topography_top_index_ji( j, i, 's' )
-               DO  k = k_wall, nzt+1
-                  f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
+
+               IF ( ijkfc(n,m,l) /= 0 )  THEN
+                  cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+                  rcorr     = fc(n,m,l) - cellsum
+                  cellsumd  = cellsumd / REAL( ijkfc(n,m,l), KIND=wp )
+               ELSE
+                  cellsum   = 0.0_wp                  
+                  rcorr     = 0.0_wp
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+! 
+!--            Distribute the correction term to the child nodes according to 
+!--            their relative difference to the parent value such that the
+!--            node with the largest difference gets the largest share of the 
+!--            correction. The distribution is skipped if rcorr is negligibly 
+!--            small in order to avoid division by zero. 
+               IF ( ABS(rcorr) < 0.000001_wp )  THEN                  
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+            
+               ijk = 1           
+               DO  i = ifl(l), ifu(l)
+                  DO  ja = jfl(m), jfu(m)
+                     DO  k = kfl(n), kfu(n)
+                        rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+                        f(k,ja,i) = f(k,ja,i) + rcorr_ijk
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
                ENDDO
-            ENDDO
-         ENDIF
-      ENDIF
+               
+            ENDDO  ! n
+         ENDDO  ! l 
+
+      ENDIF  ! var not u or w
 !
 !--   Store the boundary values also into the other redundant ghost node layers.
-!--   Please note, in case of only one ghost node layer, e.g. for the PW 
+!--   Note that in case of only one ghost node layer, e.g. for the PW 
 !--   scheme, the following loops will not be entered.
       IF ( edge == 's' )  THEN
@@ -5551,6 +5976,7 @@
  
 
-   SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z,  &
-                                  r2z, var )
+   SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,       &
+                                  r1z, r2z, kct, ifl, ifu, jfl, jfu, kfl, kfu, &
+                                  ijkfc, var )
 
 !
@@ -5576,4 +6002,16 @@
       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) ::  kc    !<
       
+      INTEGER(iwp) :: kct
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH      INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH      INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+      INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+
       CHARACTER(LEN=1), INTENT(IN) :: var   !<
 
@@ -5581,12 +6019,17 @@
       INTEGER(iwp) ::  ib  !<
       INTEGER(iwp) ::  ie  !<
+      INTEGER(iwp) ::  ijk !<
       INTEGER(iwp) ::  j   !<
-      INTEGER(iwp) ::  jb   !<
-      INTEGER(iwp) ::  je   !<      
+      INTEGER(iwp) ::  jb  !<
+      INTEGER(iwp) ::  je  !<      
       INTEGER(iwp) ::  k   !<
+      INTEGER(iwp) ::  ka  !<
       INTEGER(iwp) ::  l   !<
       INTEGER(iwp) ::  m   !<
       INTEGER(iwp) ::  n   !<
+      INTEGER(iwp) ::  var_flag  !<
       
+      REAL(wp) ::  cellsum     !<
+      REAL(wp) ::  cellsumd    !<
       REAL(wp) ::  fk          !<
       REAL(wp) ::  fkj         !<
@@ -5595,6 +6038,8 @@
       REAL(wp) ::  fkpjp       !<
       REAL(wp) ::  fkp         !<
-
-      
+      REAL(wp) ::  rcorr       !<
+      REAL(wp) ::  rcorr_ijk   !<
+
+
       IF ( var == 'w' )  THEN
          k  = nzt
@@ -5610,12 +6055,13 @@
       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
+
+      IF ( var == 'u' )  THEN 
+         var_flag = 1 
+      ELSEIF ( var == 'v' )  THEN
+         var_flag = 2 
+      ELSEIF ( var == 'w' )  THEN
+         var_flag = 3
+      ELSE
+         var_flag = 0
       ENDIF
          
@@ -5635,18 +6081,68 @@
       ENDDO
 !
-!--   Just fill up the second ghost-node layer for w.
+!--   Apply the reversibility correction to the boundary-normal velocity-
+!--   component w and scalars. It must not be applied to the boundary-
+!--   tangential velocity components u and v because their 2-D anterpolation 
+!--   cells do not cover all the child-grid nodes on the boundary.
+      IF ( .NOT. ( ( var == 'u' ) .OR. ( var == 'v' ) ) )  THEN  
+         IF ( var == 'w' )  THEN
+            n = kc(k)
+         ELSE
+            n = kc(k) + 1
+         ENDIF
+
+         DO  l = icl, icr
+            DO  m = jcs, jcn
+               ijk = 1
+               cellsum   = 0.0_wp
+               cellsumd  = 0.0_wp
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     DO  ka = kfl(n), kfu(n)
+                        cellsum = cellsum + MERGE( f(ka,j,i), 0.0_wp,           &
+                             BTEST( wall_flags_0(ka,j,i), var_flag ) )
+                        celltmpd(ijk) = ABS( fc(n,m,l) - f(ka,j,i) )
+                        cellsumd      = cellsumd  + MERGE( celltmpd(ijk),       &
+                             0.0_wp, BTEST( wall_flags_0(ka,j,i), var_flag ) )                     
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
+               ENDDO
+
+               IF ( ijkfc(n,m,l) /= 0 )  THEN
+                  cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+                  rcorr     = fc(n,m,l) - cellsum
+                  cellsumd  = cellsumd / REAL( ijkfc(n,m,l), KIND=wp )
+               ELSE
+                  cellsum   = 0.0_wp
+                  rcorr     = 0.0_wp
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+
+               IF ( ABS(rcorr) < 0.000001_wp )  THEN                  
+                  cellsumd  = 1.0_wp
+                  celltmpd  = 1.0_wp
+               ENDIF
+
+               ijk = 1
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     DO  ka = kfl(n), kfu(n)
+                        rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+                        f(ka,j,i) = f(ka,j,i) + rcorr_ijk
+                        ijk = ijk + 1
+                     ENDDO
+                  ENDDO
+               ENDDO
+               
+            ENDDO  ! m
+         ENDDO  ! l
+
+      ENDIF  ! var not u or v
+!
+!--   Just fill up the redundant second ghost-node layer for w.
       IF ( var == 'w' )  THEN
          f(nzt+1,:,:) = f(nzt,:,:)
-      ENDIF
-!
-!--   Rescale if f is the TKE.
-!--   It is assumed that the bottom surface never reaches the top boundary of a
-!--   nest domain.
-      IF ( var == 'e' )  THEN
-         DO  i = nxl, nxr
-            DO  j = nys, nyn
-               f(k,j,i) = tkefactor_t(j,i) * f(k,j,i)
-            ENDDO
-         ENDDO
       ENDIF
 
@@ -5667,15 +6163,18 @@
        CHARACTER(LEN=*), INTENT(IN) ::  var   !< identifyer for treated variable
 
+!AH       INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+       INTEGER(iwp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+
        INTEGER(iwp) ::  i         !< Running index x-direction - fine-grid
+       INTEGER(iwp) ::  icla      !< Left boundary index for anterpolation along x
+       INTEGER(iwp) ::  icra      !< Right boundary index for anterpolation along x
        INTEGER(iwp) ::  ii        !< Running index x-direction - coarse grid
-       INTEGER(iwp) ::  iclp      !< Left boundary index for anterpolation along x
-       INTEGER(iwp) ::  icrm      !< Right boundary index for anterpolation along x
        INTEGER(iwp) ::  j         !< Running index y-direction - fine-grid
+       INTEGER(iwp) ::  jcna      !< North boundary index for anterpolation along y
+       INTEGER(iwp) ::  jcsa      !< South boundary index for anterpolation along y
        INTEGER(iwp) ::  jj        !< Running index y-direction - coarse grid
-       INTEGER(iwp) ::  jcnm      !< North boundary index for anterpolation along y
-       INTEGER(iwp) ::  jcsp      !< South boundary index for anterpolation along y
        INTEGER(iwp) ::  k         !< Running index z-direction - fine-grid     
+       INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
        INTEGER(iwp) ::  kk        !< Running index z-direction - coarse grid
-       INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
        INTEGER(iwp) ::  var_flag  !< bit number used to flag topography on respective grid
 
@@ -5686,8 +6185,8 @@
        INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
        INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
-       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
-       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
-
-       INTEGER(iwp), DIMENSION(0:kct,jcs:jcn,icl:icr), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+!AH       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfl !< Indicates start index of child cells belonging to certain parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:cg%nz+1), INTENT(IN)   ::  kfu !< Indicates start index of child cells belonging to certain parent cell - z direction
 
        REAL(wp) ::  cellsum   !< sum of respective child cells belonging to parent cell 
@@ -5698,12 +6197,5 @@
   
 !
-!--    Initialize the index bounds for anterpolation 
-       iclp = icl 
-       icrm = icr 
-       jcsp = jcs 
-       jcnm = jcn 
-       kcb  = 0
-!
-!--    Define the index bounds iclp, icrm, jcsp and jcnm.
+!--    Define the index bounds icla, icra, jcsa and jcna.
 !--    Note that kcb is simply zero and kct enters here as a parameter and it is
 !--    determined in pmci_init_anterp_tophat.
@@ -5714,40 +6206,36 @@
 !--    child domain, leading to increased velocity variances. A more 
 !--    comprehensive explanation for this is still pending.
-       IF ( nesting_mode == 'vertical' )  THEN
-          IF ( bc_dirichlet_l )  THEN
-             iclp = icl + nhll
-          ENDIF
-          IF ( bc_dirichlet_r ) THEN
-             icrm = icr - nhlr
-          ENDIF
-          IF ( bc_dirichlet_s )  THEN
-             jcsp = jcs + nhls
-          ENDIF
-          IF ( bc_dirichlet_n )  THEN
-             jcnm = jcn - nhln
-          ENDIF
-       ELSE
+       icla = coarse_bound_anterp(1)
+       icra = coarse_bound_anterp(2) 
+       jcsa = coarse_bound_anterp(3)
+       jcna = coarse_bound_anterp(4)
+       kcb  = 0
+       IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
              IF ( var == 'u' )  THEN
-                iclp = icl + nhll + 1 + 1
+                icla = coarse_bound_anterp(1) + 2
              ELSE
-                iclp = icl + nhll + 1
+                icla = coarse_bound_anterp(1) + 1
              ENDIF
           ENDIF
           IF ( bc_dirichlet_r )  THEN
-             icrm = icr - nhlr - 1
+             icra = coarse_bound_anterp(2) - 1
           ENDIF
 
           IF ( bc_dirichlet_s )  THEN
              IF ( var == 'v' )  THEN
-                jcsp = jcs + nhls + 1 + 1
+                jcsa = coarse_bound_anterp(3) + 2
              ELSE
-                jcsp = jcs + nhls + 1
+                jcsa = coarse_bound_anterp(3) + 1
              ENDIF
           ENDIF
           IF ( bc_dirichlet_n )  THEN
-             jcnm = jcn - nhln - 1
+             jcna = coarse_bound_anterp(4) - 1
           ENDIF
        ENDIF
+
+!       write(9,"('pmci_anterp_tophat: ',4(e12.5,2x))")   &
+!            cg%coord_x(icla), cg%coord_y(jcsa),  cg%coord_x(icra), cg%coord_y(jcna)  
+!       flush(9)
 !
 !--    Set masking bit for topography flags 
@@ -5761,15 +6249,13 @@
           var_flag = 0
        ENDIF  
-
 !
 !--    Note that ii, jj, and kk are coarse-grid indices and i,j, and k 
 !--    are fine-grid indices.
-       DO  ii = iclp, icrm
-          DO  jj = jcsp, jcnm
+       DO  ii = icla, icra
+          DO  jj = jcsa, jcna
 !
 !--          For simplicity anterpolate within buildings and under elevated
 !--          terrain too
-             DO  kk = kcb, kct - 1
-
+             DO  kk = kcb, kct - 1                
                 cellsum = 0.0_wp
                 DO  i = ifl(ii), ifu(ii)
@@ -5783,5 +6269,7 @@
 !
 !--             Spatial under-relaxation.
-                fra  = frax(ii) * fray(jj) * fraz(kk)
+!--             The relaxation buffer zones are no longer needed with
+!--             the new reversible interpolation algorithm. 23.19.2018.
+!                fra  = frax(ii) * fray(jj) * fraz(kk)                
 !
 !--             In case all child grid points are inside topography, i.e. 
@@ -5791,8 +6279,9 @@
 !--             zero, keep the parent solution at this point. 
                 IF ( ijkfc(kk,jj,ii) /= 0 )  THEN
-                   fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +            &
-                                    fra * cellsum  /                           &
-                                    REAL( ijkfc(kk,jj,ii), KIND=wp )
-                ENDIF  
+!                   fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +         &
+!                        fra * cellsum  /                                    &
+!                        REAL( ijkfc(kk,jj,ii), KIND=wp )
+                   fc(kk,jj,ii) = cellsum / REAL( ijkfc(kk,jj,ii), KIND=wp )
+                ENDIF
 
              ENDDO
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3483)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3484)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! pmci_ensure_nest_mass_conservation is premanently removed
+! 
+! 3473 2018-10-30 20:50:15Z suehring
 ! new module for virtual measurements introduced
 ! 
@@ -517,5 +520,5 @@
     USE pmc_interface,                                                         &
         ONLY:  nested_run, nesting_mode, pmci_boundary_conds, pmci_datatrans,  &
-               pmci_ensure_nest_mass_conservation, pmci_synchronize
+               pmci_synchronize
 
     USE progress_bar,                                                          &
@@ -978,13 +981,4 @@
 !--          Set boundary conditions again after interpolation and anterpolation.
              CALL pmci_boundary_conds
-!
-!--          Correct the w top-BC in nest domains to ensure mass conservation.
-!--          This action must never be done for the root domain. Vertical 
-!--          Commented out April 18, 2018 as seemingly unnecessary. 
-!--          Will later be completely removed.
-!--             IF ( child_domain )  THEN
-!--                CALL pmci_ensure_nest_mass_conservation
-!--             ENDIF
-
 
              CALL cpu_log( log_point(60), 'nesting', 'stop' )
