Changeset 3681 for palm

Timestamp:

Jan 18, 2019 3:06:05 PM (6 years ago)

Author:

hellstea

Message:

Major update of pmc_interface_mod

Location:

palm/trunk

Files:

• SOURCE/pmc_interface_mod.f90 (modified) (133 diffs)
• UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_cbm4/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_salsagas/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 (modified) (1 diff)
• UTIL/inifor/tests/test-input-files.f90 (modified) (2 diffs)

palm/trunk/SOURCE/pmc_interface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Linear interpolations are replaced by first order interpolations. The linear 
+! interpolation routines are still included but not called. In the child 
+! inititialization the interpolation is also changed to 1st order and the linear
+! interpolation is not kept.
+! Subroutine pmci_map_fine_to_coarse_grid is rewritten.
+! Several changes in pmci_init_anterp_tophat.
+! Child's parent-grid arrays (uc, vc,...) are made non-overlapping on the PE-
+! subdomain boundaries in order to allow grid-spacing ratios higher than nbgp. 
+! Subroutine pmci_init_tkefactor is removed as unnecessary.
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Remove unused variable simulated_time
 ! 
@@ -334 +345 @@
 
 
-   USE arrays_3d,                                                              &
+    USE arrays_3d,                                                             &
         ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2,  &
                pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                   &
@@ -378 +389 @@
     USE pegrid,                                                                &
         ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
-               numprocs, pleft, pnorth, pright, psouth, status
+               numprocs, pdims, pleft, pnorth, pright, psouth, status
 
     USE pmc_child,                                                             &
@@ -446 +457 @@
 !
 !-- Geometry
-    REAL(wp), SAVE                                    ::  area_t             !<
     REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_x            !<
     REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_y            !<
@@ -455 +465 @@
 !-- 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           !<
-    REAL(wp), SAVE                              ::  xexr           !<
-    REAL(wp), SAVE                              ::  yexs           !<
-    REAL(wp), SAVE                              ::  yexn           !<
-    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !<
-    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !<
-    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !<
-    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !<
-    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !<
+    INTEGER(iwp), SAVE, DIMENSION(4), PUBLIC    ::  coarse_bound_aux    !< subdomain index bounds for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp), SAVE, DIMENSION(4), PUBLIC    ::  coarse_bound_w      !< subdomain index bounds for children's parent-grid work arrays
 
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< coarse grid array on child domain - dissipation rate
@@ -493 +494 @@
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !<
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !<
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  celltmpd !<
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !<
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !<
@@ -505 +507 @@
     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
@@ -547 +549 @@
     REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !<
 !
-!-- Upper bounds for k in anterpolation.
-    INTEGER(iwp), SAVE ::  kcto   !<
-    INTEGER(iwp), SAVE ::  kctw   !<
+    INTEGER(iwp), SAVE ::  igsr     !< Integer grid-spacing ratio in i-direction 
+    INTEGER(iwp), SAVE ::  jgsr     !< Integer grid-spacing ratio in j-direction 
+    INTEGER(iwp), SAVE ::  kgsr     !< Integer grid-spacing ratio in k-direction 
+    INTEGER(iwp), SAVE ::  kcto     !< Upper bound for k in anterpolation of variables other than w.
+    INTEGER(iwp), SAVE ::  kctw     !< Upper bound for k in anterpolation of w.
+    INTEGER(iwp), SAVE ::  nxlfc    !< Lower index limit in x-direction for fine-to-coarse index mapping and interpolaton coefficient arrays
+    INTEGER(iwp), SAVE ::  nxrfc    !< Upper index limit in x-direction for fine-to-coarse index mapping and interpolaton coefficient arrays
+    INTEGER(iwp), SAVE ::  nynfc    !< Upper index limit in y-direction for fine-to-coarse index mapping and interpolaton coefficient arrays
+    INTEGER(iwp), SAVE ::  nysfc    !< Lower index limit in y-direction for fine-to-coarse index mapping and interpolaton coefficient arrays
 !
 !-- Upper bound for k in log-law correction in interpolation.
@@ -556 +564 @@
     INTEGER(iwp), SAVE ::  nzt_topo_nestbc_r   !<
     INTEGER(iwp), SAVE ::  nzt_topo_nestbc_s   !<
-!
-!-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
-    INTEGER(iwp), SAVE ::  nhll   !<
-    INTEGER(iwp), SAVE ::  nhlr   !<
-    INTEGER(iwp), SAVE ::  nhls   !<
-    INTEGER(iwp), SAVE ::  nhln   !<
 !
 !-- Spatial under-relaxation coefficients for anterpolation.
@@ -633 +635 @@
     INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: nr_part  !<
     INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: part_adr !<
+
+!AH 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_lr
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_sn
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_t
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarrc_lr
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarrc_sn
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarrc_t
+    INTEGER(iwp) :: workarrc_lr_exchange_type
+    INTEGER(iwp) :: workarrc_sn_exchange_type
+    INTEGER(iwp) :: workarrc_t_exchange_type_x
+    INTEGER(iwp) :: workarrc_t_exchange_type_y
+!AH
     
     INTERFACE pmci_boundary_conds
@@ -914 +929 @@
           nz_cl = nz
 !
-!--       Find the highest nest level in the coarse grid for the reduced z
+!--       Find the highest nest level in the parent grid for the reduced z
 !--       transfer
           DO  k = 1, nz                 
@@ -1130 +1145 @@
           CALL MPI_COMM_SIZE( comm1dy, npy, ierr )
 !
-!--       The +1 in index is because PALM starts with nx=0
+!--       Nrx is the same for all PEs and so is nry, thus there is no need to compute
+!--       them separately for each PE.
           nrx = nxr - nxl + 1
           nry = nyn - nys + 1
@@ -1139 +1155 @@
 !
 !--          Area along y required by actual child PE
-             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
+             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)  !: j = jcs, jcn of PE# k 
 !
 !--             Area along x required by actual child PE
-                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
+                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)  !: i = icl, icr of PE# k 
 
                    px = i / nrx
@@ -1230 +1246 @@
     INTEGER(iwp) ::  i          !<
     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) ::  icl        !< Left index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  icla       !< Left index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp) ::  iclw       !< Left index limit for children's parent-grid work arrays
+    INTEGER(iwp) ::  icr        !< Left index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  icra       !< Right index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp) ::  icrw       !< Right index limit for children's parent-grid work arrays
     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
+    INTEGER(iwp) ::  jcn        !< North index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  jcna       !< North index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp) ::  jcnw       !< North index limit for children's parent-grid work arrays
+    INTEGER(iwp) ::  jcs        !< South index limit for children's parent-grid arrays
+    INTEGER(iwp) ::  jcsa       !< South index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp) ::  jcsw       !< South index limit for children's parent-grid work arrays
+    INTEGER(iwp) ::  n          !< Running index for number of chemical species
 
     INTEGER(iwp), DIMENSION(5) ::  val        !<
@@ -1441 +1461 @@
        ENDIF
 !
-!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio. Only
-!--    if both parent and child are in LES mode or in RANS mode. 
-!--    Please note, in case parent and child are in RANS mode, TKE weighting 
-!--    factor is simply one. 
-       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
-            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
-               .NOT. constant_diffusion ) )  CALL pmci_init_tkefactor
-!
 !--    Two-way coupling for general and vertical nesting.
 !--    Precompute the index arrays and relaxation functions for the
@@ -1457 +1469 @@
 !--    included in the interpolation algorithms.
        CALL pmci_init_anterp_tophat
-!
-!--    Finally, compute the total area of the top-boundary face of the domain.
-!--    This is needed in the pmc_ensure_nest_mass_conservation      
-       area_t = ( nx + 1 ) * (ny + 1 ) * dx * dy
 
     ENDIF
@@ -1476 +1484 @@
        INTEGER(iwp), DIMENSION(2)          ::  size_of_array      !<
                                    
-       INTEGER(iwp) :: i        !<   
+       INTEGER(iwp) :: i        !< 
+       INTEGER(iwp) :: iauxl    !<
+       INTEGER(iwp) :: iauxr    !<
        INTEGER(iwp) :: ijaux    !< 
-       INTEGER(iwp) :: j        !< 
+       INTEGER(iwp) :: j        !<
+       INTEGER(iwp) :: jauxs    !<
+       INTEGER(iwp) :: jauxn    !< 
        REAL(wp) ::  loffset     !<
        REAL(wp) ::  noffset     !<
        REAL(wp) ::  roffset     !<
        REAL(wp) ::  soffset     !<
-
-!
-!--    If the fine- and coarse grid nodes do not match:
-       loffset = MOD( coord_x(nxl), cg%dx )
-       xexl    = cg%dx + loffset
-!
-!--    This is needed in the anterpolation phase
-       nhll = CEILING( xexl / cg%dx )
-       xcs  = coord_x(nxl) - xexl
-       DO  i = 0, cg%nx
-          IF ( cg%coord_x(i) > xcs )  THEN
-             icl = MAX( -1, i-1 )
-             EXIT
-          ENDIF
-       ENDDO
-!
-!--    If the fine- and coarse grid nodes do not match
-       roffset = MOD( coord_x(nxr+1), cg%dx )
-       xexr    = cg%dx + roffset
-!
-!--    This is needed in the anterpolation phase
-       nhlr = CEILING( xexr / cg%dx )
-       xce  = coord_x(nxr+1) + xexr
-!--    One "extra" layer is taken behind the right boundary 
-!--    because it may be needed in cases of non-integer grid-spacing ratio
-       DO  i = cg%nx, 0 , -1
-          IF ( cg%coord_x(i) < xce )  THEN
-             icr = MIN( cg%nx+1, i+1 )
-             EXIT
-          ENDIF
-       ENDDO
-!
-!--    If the fine- and coarse grid nodes do not match
-       soffset = MOD( coord_y(nys), cg%dy )
-       yexs    = cg%dy + soffset
-!
-!--    This is needed in the anterpolation phase
-       nhls = CEILING( yexs / cg%dy )
-       ycs  = coord_y(nys) - yexs
-       DO  j = 0, cg%ny
-          IF ( cg%coord_y(j) > ycs )  THEN
-             jcs = MAX( -nbgp, j-1 )
-             EXIT
-          ENDIF
-       ENDDO
-!
-!--    If the fine- and coarse grid nodes do not match
-       noffset = MOD( coord_y(nyn+1), cg%dy )
-       yexn    = cg%dy + noffset
-!
-!--    This is needed in the anterpolation phase
-       nhln = CEILING( yexn / cg%dy )
-       yce  = coord_y(nyn+1) + yexn
-!--    One "extra" layer is taken behind the north boundary 
-!--    because it may be needed in cases of non-integer grid-spacing ratio
-       DO  j = cg%ny, 0, -1
-          IF ( cg%coord_y(j) < yce )  THEN
-             jcn = MIN( cg%ny + nbgp, j+1 )
-             EXIT
-          ENDIF
-       ENDDO
-
-       coarse_bound(1) = icl
-       coarse_bound(2) = icr
-       coarse_bound(3) = jcs
-       coarse_bound(4) = jcn
-       coarse_bound(5) = myid
+       REAL(wp) ::  xexl        !< Parent-grid array exceedance behind the left edge of the child PE subdomain
+       REAL(wp) ::  xexr        !< Parent-grid array exceedance behind the right edge of the child PE subdomain 
+       REAL(wp) ::  yexs        !< Parent-grid array exceedance behind the south edge of the child PE subdomain
+       REAL(wp) ::  yexn        !< Parent-grid array exceedance behind the north edge of the child PE subdomain
+
+!AH!
+!AH!--    If the fine- and coarse grid nodes do not match:
+!AH       loffset = MOD( coord_x(nxl), cg%dx )
+!AH       xexl    = cg%dx + loffset
+!AH       xcs  = coord_x(nxl) - xexl
+!AH       DO  i = 0, cg%nx
+!AH          IF ( cg%coord_x(i) > xcs )  THEN
+!AH             icl = MAX( -1, i-1 )
+!AH             EXIT
+!AH          ENDIF
+!AH       ENDDO
+!AH!
+!AH!--    If the fine- and coarse grid nodes do not match
+!AH       roffset = MOD( coord_x(nxr+1), cg%dx )
+!AH       xexr    = cg%dx + roffset
+!AH       xce  = coord_x(nxr+1)
+!AH       IF ( nxr == nx )  THEN
+!AH          xce = xce + xexr
+!AH       ENDIF
+!AH       DO  i = cg%nx, 0 , -1
+!AH          IF ( cg%coord_x(i) < xce )  THEN
+!AH             icr = MIN( cg%nx+1, i+1 )
+!AH             EXIT
+!AH          ENDIF
+!AH       ENDDO
+!AH!
+!AH!--    If the fine- and coarse grid nodes do not match
+!AH       soffset = MOD( coord_y(nys), cg%dy )
+!AH       yexs    = cg%dy + soffset
+!AH       ycs  = coord_y(nys) - yexs
+!AH       DO  j = 0, cg%ny
+!AH          IF ( cg%coord_y(j) > ycs )  THEN
+!AH             jcs = MAX( -nbgp, j-1 )
+!AH             EXIT
+!AH          ENDIF
+!AH       ENDDO
+!AH!
+!AH!--    If the fine- and coarse grid nodes do not match
+!AH       noffset = MOD( coord_y(nyn+1), cg%dy )
+!AH       yexn    = cg%dy + noffset
+!AH       yce  = coord_y(nyn+1)
+!AH       IF ( nyn == ny )  THEN
+!AH          yce = yce + yexn
+!AH       ENDIF
+!AH       DO  j = cg%ny, 0, -1
+!AH          IF ( cg%coord_y(j) < yce )  THEN
+!AH             jcn = MIN( cg%ny + nbgp, j+1 )
+!AH             EXIT
+!AH          ENDIF
+!AH       ENDDO
+!AH
+!AH       coarse_bound(1) = icl
+!AH       coarse_bound(2) = icr
+!AH       coarse_bound(3) = jcs
+!AH       coarse_bound(4) = jcn
+!AH       coarse_bound(5) = myid
 !
 !--    Determine the anterpolation index limits. If at least half of the 
@@ -1557 +1563 @@
 !--    anterpolation domain, or not included at all if we are at the outer 
 !--    edge of the child domain.
+
+!
+!--    Left
+       IF  ( bc_dirichlet_l )  THEN
+          loffset = MOD( coord_x(nxl), cg%dx )
+          xexl  = 2 * cg%dx + loffset
+          iauxl = 0
+       ELSE
+          xexl  = 0.0_wp
+          iauxl = 1
+       ENDIF
+       xcs     = coord_x(nxl) - xexl
        DO  i = 0, cg%nx
-          IF ( cg%coord_x(i) + 0.5_wp * cg%dx >= coord_x(nxl) )  THEN
-             icla = MAX( 0, i )
+          IF ( cg%coord_x(i) + 0.5_wp * cg%dx >= xcs )  THEN
+             icl = MAX( 0, i )
              EXIT
           ENDIF
        ENDDO
+!
+!--    Right
+       IF  ( bc_dirichlet_r )  THEN
+          roffset = MOD( coord_x(nxr+1), cg%dx )
+          xexr  = 2 * cg%dx + roffset
+          iauxr = 0  
+       ELSE
+          xexr  = 0.0_wp
+          iauxr = 1  
+       ENDIF
+       xce  = coord_x(nxr+1) + xexr
        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 )
+          IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= xce )  THEN
+             icr = MIN( cg%nx, i )
              EXIT
           ENDIF
        ENDDO
+!
+!--    South
+       IF  ( bc_dirichlet_s )  THEN
+          soffset = MOD( coord_y(nys), cg%dy )
+          yexs  = 2 * cg%dy + soffset
+          jauxs = 0  
+       ELSE
+          yexs  = 0.0_wp
+          jauxs = 1  
+       ENDIF
+       ycs  = coord_y(nys) - yexs
        DO  j = 0, cg%ny
-          IF ( cg%coord_y(j) + 0.5_wp * cg%dy >= coord_y(nys) )  THEN
-             jcsa = MAX( 0, j )
+          IF ( cg%coord_y(j) + 0.5_wp * cg%dy >= ycs )  THEN
+             jcs = MAX( 0, j )
              EXIT
           ENDIF
        ENDDO
+!
+!--    North
+       IF  ( bc_dirichlet_n )  THEN
+          noffset = MOD( coord_y(nyn+1), cg%dy )
+          yexn  = 2 * cg%dy + noffset
+          jauxn = 0
+       ELSE
+          yexn  = 0.0_wp
+          jauxn = 1
+       ENDIF
+       yce  = coord_y(nyn+1) + yexn
        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 )
+          IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= yce )  THEN
+             jcn = MIN( cg%ny, j )
              EXIT
           ENDIF
        ENDDO
 !
-!--    Make sure that the indexing is contiguous 
+!--    Make sure that the indexing is contiguous (no gaps, no overlaps) 
+#if defined( __parallel )
        IF ( nxl == 0 )  THEN
-          CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+          CALL MPI_SEND( icr, 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
+          icl = ijaux + 1
        ELSE
-          CALL MPI_SEND( icra, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+          CALL MPI_SEND( icr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 
-          icla = ijaux + 1
+          icl = ijaux + 1
        ENDIF
        IF ( nys == 0 )  THEN
-          CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+          CALL MPI_SEND( jcn, 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
+          jcs = ijaux + 1
        ELSE
-          CALL MPI_SEND( jcna, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+          CALL MPI_SEND( jcn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 
-          jcsa = ijaux + 1
+          jcs = 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
+#endif       
+
+       WRITE(9,"('Pmci_map_fine_to_coarse_grid. parent-grid array bounds: ',4(i3,2x))") icl, icr, jcs, jcn
+       FLUSH(9)
+
+       coarse_bound(1) = icl
+       coarse_bound(2) = icr
+       coarse_bound(3) = jcs
+       coarse_bound(4) = jcn
+       coarse_bound(5) = myid
+!
+!--    The following index bounds are used for allocating index mapping and some other auxiliary arrays
+       coarse_bound_aux(1) = icl - iauxl
+       coarse_bound_aux(2) = icr + iauxr
+       coarse_bound_aux(3) = jcs - jauxs
+       coarse_bound_aux(4) = jcn + jauxn
 !
 !--    Note that MPI_Gather receives data from all processes in the rank order
@@ -1636 +1697 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  acsize  !<
-       INTEGER(iwp) ::  i       !<
-       INTEGER(iwp) ::  j       !<
-       INTEGER(iwp) ::  k       !<
-       INTEGER(iwp) ::  kc      !<
-       INTEGER(iwp) ::  kdzo    !<
-       INTEGER(iwp) ::  kdzw    !<       
-
-       REAL(wp) ::  dzmin       !<
-       REAL(wp) ::  parentdzmax !<
-       REAL(wp) ::  xb          !<
-       REAL(wp) ::  xcsu        !<
-       REAL(wp) ::  xfso        !<
-       REAL(wp) ::  xcso        !<
-       REAL(wp) ::  xfsu        !<
-       REAL(wp) ::  yb          !<
-       REAL(wp) ::  ycso        !<
-       REAL(wp) ::  ycsv        !<
-       REAL(wp) ::  yfso        !<
-       REAL(wp) ::  yfsv        !<
-       REAL(wp) ::  zcso        !<
-       REAL(wp) ::  zcsw        !<
-       REAL(wp) ::  zfso        !<
-       REAL(wp) ::  zfsw        !<
+       INTEGER(iwp) ::  acsize   !< Maximum dimension of anterpolation cell.
+       INTEGER(iwp) ::  i        !< Child-grid i-index
+       INTEGER(iwp) ::  ierr     !< MPI error code
+       INTEGER(iwp) ::  j        !< Child-grid j-index
+       INTEGER(iwp) ::  k        !< Child-grid k-index
+       INTEGER(iwp) ::  kc       !< 
+       INTEGER(iwp) ::  kdzo     !<
+       INTEGER(iwp) ::  kdzw     !<
+       INTEGER(iwp) ::  moff     !< Parent-grid bound offset in j-direction
+       INTEGER(iwp) ::  loff     !< Parent-grid bound offset in i-direction
+
+       REAL(wp) ::  dzmin        !<
+       REAL(wp) ::  parentdzmax  !<
+       REAL(wp) ::  xb           !<
+       REAL(wp) ::  xcsu         !<
+       REAL(wp) ::  xfso         !<
+       REAL(wp) ::  xcso         !<
+       REAL(wp) ::  xfsu         !<
+       REAL(wp) ::  yb           !<
+       REAL(wp) ::  ycso         !<
+       REAL(wp) ::  ycsv         !<
+       REAL(wp) ::  yfso         !<
+       REAL(wp) ::  yfsv         !<
+       REAL(wp) ::  zcso         !<
+       REAL(wp) ::  zcsw         !<
+       REAL(wp) ::  zfso         !<
+       REAL(wp) ::  zfsw         !<
       
        
-       xb = nxl * dx
-       yb = nys * dy
-      
-       ALLOCATE( icu(nxlg:nxrg) )
-       ALLOCATE( ico(nxlg:nxrg) )
-       ALLOCATE( jcv(nysg:nyng) )
-       ALLOCATE( jco(nysg:nyng) )
-       ALLOCATE( kcw(nzb:nzt+1) )
-       ALLOCATE( kco(nzb:nzt+1) )
-       ALLOCATE( r1xu(nxlg:nxrg) )
-       ALLOCATE( r2xu(nxlg:nxrg) )
-       ALLOCATE( r1xo(nxlg:nxrg) )
-       ALLOCATE( r2xo(nxlg:nxrg) )
-       ALLOCATE( r1yv(nysg:nyng) )
-       ALLOCATE( r2yv(nysg:nyng) )
-       ALLOCATE( r1yo(nysg:nyng) )
-       ALLOCATE( r2yo(nysg:nyng) )
-       ALLOCATE( r1zw(nzb:nzt+1) )
-       ALLOCATE( r2zw(nzb:nzt+1) )
-       ALLOCATE( r1zo(nzb:nzt+1) )
-       ALLOCATE( r2zo(nzb:nzt+1) )
-!
-!--    Note that the node coordinates xfs... and xcs... are relative to the
-!--    lower-left-bottom corner of the fc-array, not the actual child domain
-!--    corner
-       DO  i = nxlg, nxrg
-          xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
-          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
-          icu(i)  = icl + FLOOR( xfsu / cg%dx )
-          ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * cg%dx ) / cg%dx )
-          xcsu    = ( icu(i) - icl ) * cg%dx
-          xcso    = ( ico(i) - icl ) * cg%dx + 0.5_wp * cg%dx
+!AH
+!
+!--    Allocate child-grid work arrays for interpolation.
+       CALL pmci_allocate_finegrid_workarrays
+!       
+!--    Determine index bounds for the parent-grid work arrays for
+!--    interpolation and allocate them.
+       CALL pmci_allocate_coarsegrid_workarrays
+!       
+!--    Define the MPI-datatypes for parent-grid work array 
+!--    exchange between the PE-subdomains.
+       CALL pmci_create_coarsegrid_workarray_exchange_datatypes
+!
+!--    Determine index bounds for the fine-to-coarse grid index mapping arrays 
+!--    and interpolation-coefficient arrays and allocate them.       
+       CALL pmci_allocate_fine_to_coarse_mapping_arrays
+!AH       
+!
+       xb   = nxl * dx
+       IF  ( bc_dirichlet_l )  THEN
+          loff = 2
+       ELSE
+          loff = 0
+       ENDIF
+!AH       DO  i = nxlg, nxrg
+       DO  i = nxl-1, nxr+1
+          xfsu    = coord_x(i) - ( lower_left_coord_x + xb )
+          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb )
+!
+!--       icl points to 2 parent-grid cells left form the left nest boundary,
+!--       thence icl + loff points to the left nest boundary.
+          icu(i)  = icl + loff + FLOOR( xfsu / cg%dx )
+          ico(i)  = icl + loff + FLOOR( ( xfso - 0.5_wp * cg%dx ) / cg%dx )
+          xcsu    = ( icu(i) - ( icl + loff ) ) * cg%dx
+          xcso    = ( ico(i) - ( icl + loff ) + 0.5_wp ) * cg%dx
           r2xu(i) = ( xfsu - xcsu ) / cg%dx
           r2xo(i) = ( xfso - xcso ) / cg%dx
@@ -1699 +1766 @@
           r1xo(i) = 1.0_wp - r2xo(i)
        ENDDO
-
-       DO  j = nysg, nyng
-          yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
-          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
-          jcv(j)  = jcs + FLOOR( yfsv / cg%dy )
-          jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * cg%dy ) / cg%dy )
-          ycsv    = ( jcv(j) - jcs ) * cg%dy
-          ycso    = ( jco(j) - jcs ) * cg%dy + 0.5_wp * cg%dy
+!       
+!--    Fill up the values behind nest boundaries by copying from inside 
+!--    the domain.
+       IF  ( bc_dirichlet_l )  THEN
+          icu(nxlfc:nxl-1)  = icu(nxlfc+igsr:nxl-1+igsr) - 1
+          r1xu(nxlfc:nxl-1) = r1xu(nxlfc+igsr:nxl-1+igsr)
+          r2xu(nxlfc:nxl-1) = r2xu(nxlfc+igsr:nxl-1+igsr)
+          ico(nxlfc:nxl-1)  = ico(nxlfc+igsr:nxl-1+igsr) - 1
+          r1xo(nxlfc:nxl-1) = r1xo(nxlfc+igsr:nxl-1+igsr)
+          r2xo(nxlfc:nxl-1) = r2xo(nxlfc+igsr:nxl-1+igsr)
+       ENDIF
+       
+       IF  ( bc_dirichlet_r )  THEN
+          icu(nxr+1:nxrfc)  = icu(nxr+1-igsr:nxrfc-igsr) + 1
+          r1xu(nxr+1:nxrfc) = r1xu(nxr+1-igsr:nxrfc-igsr)
+          r2xu(nxr+1:nxrfc) = r2xu(nxr+1-igsr:nxrfc-igsr)
+          ico(nxr+1:nxrfc)  = ico(nxr+1-igsr:nxrfc-igsr) + 1
+          r1xo(nxr+1:nxrfc) = r1xo(nxr+1-igsr:nxrfc-igsr)
+          r2xo(nxr+1:nxrfc) = r2xo(nxr+1-igsr:nxrfc-igsr)
+       ENDIF
+!
+!--    Print out the indices and coefficients for checking and debugging purposes
+       DO  i = nxlfc, nxrfc
+          WRITE(9,"('pmci_init_interp_tril: i, icu, r1xu r2xu ', 2(i4,2x),2(e12.5,2x))") &
+               i, icu(i), r1xu(i), r2xu(i)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
+       DO  i = nxlfc, nxrfc
+          WRITE(9,"('pmci_init_interp_tril: i, ico, r1xo r2xo ', 2(i4,2x),2(e12.5,2x))") &
+               i, ico(i), r1xo(i), r2xo(i)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
+
+       yb   = nys * dy
+       IF  ( bc_dirichlet_s )  THEN
+          moff = 2
+       ELSE
+          moff = 0
+       ENDIF
+!AH       DO  j = nysg, nyng
+       DO  j = nys-1, nyn+1
+          yfsv    = coord_y(j) - ( lower_left_coord_y + yb )
+          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb )
+!
+!--       jcs points to 2 parent-grid cells south form the south nest boundary,
+!--       thence jcs + moff points to the south nest boundary.
+          jcv(j)  = jcs + moff + FLOOR( yfsv / cg%dy )
+          jco(j)  = jcs + moff + FLOOR( ( yfso - 0.5_wp * cg%dy ) / cg%dy )
+          ycsv    = ( jcv(j) - ( jcs + moff ) ) * cg%dy
+          ycso    = ( jco(j) - ( jcs + moff ) + 0.5_wp ) * cg%dy
           r2yv(j) = ( yfsv - ycsv ) / cg%dy
           r2yo(j) = ( yfso - ycso ) / cg%dy
@@ -1712 +1823 @@
           r1yo(j) = 1.0_wp - r2yo(j)
        ENDDO
+!       
+!--    Fill up the values behind nest boundaries by copying from inside 
+!--    the domain.
+       IF  ( bc_dirichlet_s )  THEN
+          jcv(nysfc:nys-1)  = jcv(nysfc+jgsr:nys-1+jgsr) - 1
+          r1yv(nysfc:nys-1) = r1yv(nysfc+jgsr:nys-1+jgsr)
+          r2yv(nysfc:nys-1) = r2yv(nysfc+jgsr:nys-1+jgsr)
+          jco(nysfc:nys-1)  = jco(nysfc+jgsr:nys-1+jgsr) - 1
+          r1yo(nysfc:nys-1) = r1yo(nysfc+jgsr:nys-1+jgsr)
+          r2yo(nysfc:nys-1) = r2yo(nysfc+jgsr:nys-1+jgsr)
+       ENDIF
+       
+       IF  ( bc_dirichlet_n )  THEN
+          jcv(nyn+1:nynfc)  = jcv(nyn+1-jgsr:nynfc-jgsr) + 1
+          r1yv(nyn+1:nynfc) = r1yv(nyn+1-jgsr:nynfc-jgsr)
+          r2yv(nyn+1:nynfc) = r2yv(nyn+1-jgsr:nynfc-jgsr)
+          jco(nyn+1:nynfc)  = jco(nyn+1-jgsr:nynfc-jgsr) + 1
+          r1yo(nyn+1:nynfc) = r1yo(nyn+1-jgsr:nynfc-jgsr)
+          r2yo(nyn+1:nynfc) = r2yo(nyn+1-jgsr:nynfc-jgsr)
+       ENDIF
+!
+!--    Print out the indices and coefficients for checking and debugging purposes
+       DO  j = nysfc, nynfc
+          WRITE(9,"('pmci_init_interp_tril: j, jcv, r1yv r2yv ', 2(i4,2x),2(e12.5,2x))") &
+               j, jcv(j), r1yv(j), r2yv(j)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
+       DO  j = nysfc, nynfc
+          WRITE(9,"('pmci_init_interp_tril: j, jco, r1yo r2yo ', 2(i4,2x),2(e12.5,2x))") &
+               j, jco(j), r1yo(j), r2yo(j)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
 
        DO  k = nzb, nzt + 1
@@ -1736 +1881 @@
           r1zo(k) = 1.0_wp - r2zo(k)
        ENDDO
+!
+!--    Set the interpolation index- and coefficient-information to the 
+!--    child-grid cells within the uppermost parent-grid cell. This 
+!--    information is only needed for the reversibility correction.   
+       kco(nzt+2:nzt+kgsr)  = kco(nzt+1) + 1 
+       r1zo(nzt+2:nzt+kgsr) = r1zo(nzt+2-kgsr:nzt)
+       r2zo(nzt+2:nzt+kgsr) = r2zo(nzt+2-kgsr:nzt)
+! 
+!--    kcw, r1zw and r2zw are not needed when k > nzt+1
+       kcw(nzt+2:nzt+kgsr)  = 0
+       r1zw(nzt+2:nzt+kgsr) = 0.0_wp
+       r2zw(nzt+2:nzt+kgsr) = 0.0_wp
+!
+!--    Print out the indices and coefficients for checking and debugging purposes
+       DO  k = nzb, nzt+1
+          WRITE(9,"('pmci_init_interp_tril: k, kcw, r1zw r2zw ', 2(i4,2x),2(e12.5,2x))") &
+               k, kcw(k), r1zw(k), r2zw(k)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
+       DO  k = nzb, nzt + kgsr
+          WRITE(9,"('pmci_init_interp_tril: k, kco, r1zo r2zo ', 2(i4,2x),2(e12.5,2x))") &
+               k, kco(k), r1zo(k), r2zo(k)
+          FLUSH(9)
+       ENDDO
+       WRITE(9,*)
 !
 !--    Determine the maximum dimension of anterpolation cells and allocate the
@@ -1751 +1922 @@
             CEILING( parentdzmax / dzmin )
        ALLOCATE( celltmpd(1:acsize) )
-!       write(9,"('acsize: ',i3,2(e12.5,2x))") acsize, dzmin, parentdzmax
       
     END SUBROUTINE pmci_init_interp_tril
-
-
-
+    
+
+
+    SUBROUTINE pmci_allocate_finegrid_workarrays
+!
+!--    Allocate child-grid work-arrays for interpolation
+       IMPLICIT NONE
+
+
+       igsr = NINT( cg%dx / dx, iwp )
+       jgsr = NINT( cg%dy / dy, iwp )
+       kgsr = NINT( cg%dzu(cg%nz+1) / dzu(nzt+1), iwp )
+       write(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
+       flush(9)
+!       
+!--    Note that i-indexing for workarr_lr runs from 0 to igsr-1.
+!--    For u, only 0-element is used and all elements 0:igsr-1
+!--    are used for all other variables.
+       ALLOCATE( workarr_lr(nzb:nzt+1,nysg:nyng,0:igsr-1) )
+!       
+!--    Note that j-indexing for workarr_sn runs from 0 to jgsr-1.
+!--    For v, only 0-element is used and all elements 0:jgsr-1
+!--    are used for all other variables.
+       ALLOCATE( workarr_sn(nzb:nzt+1,0:jgsr-1,nxlg:nxrg) )
+!
+!--    Note that genuine k-indexing is used for workarr_t.
+!--    Only nzt-element is used for w and elements nzt+1:nzt+kgsr
+!--    are used for all other variables. 
+       ALLOCATE( workarr_t(nzt:nzt+kgsr,nysg:nyng,nxlg:nxrg) )
+
+    END SUBROUTINE pmci_allocate_finegrid_workarrays
+
+
+
+    SUBROUTINE pmci_allocate_coarsegrid_workarrays
+!
+!--    Allocate parent-grid work-arrays for interpolation
+       IMPLICIT NONE
+
+!
+!--    Determine and store the PE-subdomain dependent index bounds
+       IF  ( bc_dirichlet_l )  THEN
+          iclw = icl + 1
+       ELSE
+          iclw = icl - 1
+       ENDIF
+
+       IF  ( bc_dirichlet_r )  THEN
+          icrw = icr - 1
+       ELSE
+          icrw = icr + 1
+       ENDIF
+
+       IF  ( bc_dirichlet_s )  THEN
+          jcsw = jcs + 1
+       ELSE
+          jcsw = jcs - 1
+       ENDIF
+
+       IF  ( bc_dirichlet_n )  THEN
+          jcnw = jcn - 1
+       ELSE
+          jcnw = jcn + 1
+       ENDIF
+    
+       coarse_bound_w(1) = iclw
+       coarse_bound_w(2) = icrw
+       coarse_bound_w(3) = jcsw
+       coarse_bound_w(4) = jcnw
+!
+!--    Left and right boundaries.
+       ALLOCATE( workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2) )
+!
+!--    South and north boundaries.
+       ALLOCATE( workarrc_sn(0:cg%nz+1,0:2,iclw:icrw) )
+!
+!--    Top boundary.
+       ALLOCATE( workarrc_t(0:2,jcsw:jcnw,iclw:icrw) )
+
+    END SUBROUTINE pmci_allocate_coarsegrid_workarrays
+
+
+
+    SUBROUTINE pmci_create_coarsegrid_workarray_exchange_datatypes
+!
+!--    Define specific MPI types for workarrc-exhchange.
+       IMPLICIT NONE
+
+#if defined( __parallel )       
+!
+!--    For the left and right boundaries
+       CALL MPI_TYPE_VECTOR( 3, cg%nz+2, (jcnw-jcsw+1)*(cg%nz+2), MPI_REAL,     &
+            workarrc_lr_exchange_type, ierr )
+       CALL MPI_TYPE_COMMIT( workarrc_lr_exchange_type, ierr )
+!
+!--    For the south and north boundaries
+       CALL MPI_TYPE_VECTOR( 1, 3*(cg%nz+2), 3*(cg%nz+2), MPI_REAL,             &
+            workarrc_sn_exchange_type, ierr )
+       CALL MPI_TYPE_COMMIT( workarrc_sn_exchange_type, ierr )
+!
+!--    For the top-boundary x-slices
+       CALL MPI_TYPE_VECTOR( icrw-iclw+1, 3, 3*(jcnw-jcsw+1), MPI_REAL,         &
+            workarrc_t_exchange_type_x, ierr )
+       CALL MPI_TYPE_COMMIT( workarrc_t_exchange_type_x, ierr )
+!
+!--    For the top-boundary y-slices
+       CALL MPI_TYPE_VECTOR( 1, 3*(jcnw-jcsw+1), 3*(jcnw-jcsw+1), MPI_REAL,         &
+            workarrc_t_exchange_type_y, ierr )
+       CALL MPI_TYPE_COMMIT( workarrc_t_exchange_type_y, ierr )
+#endif
+       
+    END SUBROUTINE pmci_create_coarsegrid_workarray_exchange_datatypes
+
+
+
+    SUBROUTINE pmci_allocate_fine_to_coarse_mapping_arrays
+!
+!--    Define index limits and allocate the fine-to-coarse grid index mapping 
+!--    arrays and interpolation coefficient arrays.      
+       IMPLICIT NONE
+
+       
+       IF  ( bc_dirichlet_l )  THEN
+!AH          nxlfc = MIN( nxl-igsr, nxlg )          
+          nxlfc = nxl - igsr
+       ELSE
+!AH          nxlfc = nxlg
+          nxlfc = nxl - 1
+       ENDIF
+       IF  ( bc_dirichlet_r )  THEN
+!AH          nxrfc = MAX( nxr+igsr, nxrg )
+          nxrfc = nxr + igsr
+       ELSE
+!AH          nxrfc = nxrg
+          nxrfc = nxr + 1          
+       ENDIF
+
+       IF  ( bc_dirichlet_s )  THEN
+!AH          nysfc = MIN( nys-jgsr, nysg )
+          nysfc = nys - jgsr
+       ELSE
+!AH          nysfc = nysg
+          nysfc = nys - 1
+       ENDIF
+       IF  ( bc_dirichlet_n )  THEN
+!AH          nynfc = MAX( nyn+jgsr, nyng )
+          nynfc = nyn + jgsr
+       ELSE
+!AH          nynfc = nyng          
+          nynfc = nyn + 1          
+       ENDIF
+
+       ALLOCATE( icu(nxlfc:nxrfc) )       
+       ALLOCATE( ico(nxlfc:nxrfc) )
+       ALLOCATE( jcv(nysfc:nynfc) )
+       ALLOCATE( jco(nysfc:nynfc) )
+       ALLOCATE( kcw(nzb:nzt+kgsr) )
+       ALLOCATE( kco(nzb:nzt+kgsr) )
+
+       ALLOCATE( r1xu(nxlfc:nxrfc) )
+       ALLOCATE( r2xu(nxlfc:nxrfc) )
+       ALLOCATE( r1xo(nxlfc:nxrfc) )
+       ALLOCATE( r2xo(nxlfc:nxrfc) )
+       ALLOCATE( r1yv(nysfc:nynfc) )
+       ALLOCATE( r2yv(nysfc:nynfc) )
+       ALLOCATE( r1yo(nysfc:nynfc) )
+       ALLOCATE( r2yo(nysfc:nynfc) )
+       ALLOCATE( r1zw(nzb:nzt+kgsr) )
+       ALLOCATE( r2zw(nzb:nzt+kgsr) )
+       ALLOCATE( r1zo(nzb:nzt+kgsr) )
+       ALLOCATE( r2zo(nzb:nzt+kgsr) )
+
+    END SUBROUTINE pmci_allocate_fine_to_coarse_mapping_arrays
+
+     
+       
     SUBROUTINE pmci_init_loglaw_correction
 !
@@ -3092 +3435 @@
        INTEGER(iwp) ::  istart   !<
        INTEGER(iwp) ::  ir       !<
+       INTEGER(iwp) ::  iw       !< Fine-grid index limited to -1 <= iw <= nx+1
        INTEGER(iwp) ::  j        !< Fine-grid index
        INTEGER(iwp) ::  jj       !< Coarse-grid index
        INTEGER(iwp) ::  jstart   !<
        INTEGER(iwp) ::  jr       !<
+       INTEGER(iwp) ::  jw       !< Fine-grid index limited to -1 <= jw <= ny+1
        INTEGER(iwp) ::  k        !< Fine-grid index
        INTEGER(iwp) ::  kk       !< Coarse-grid index
        INTEGER(iwp) ::  kstart   !<
+       INTEGER(iwp) ::  kw       !< Fine-grid index limited to kw <= nzt+1
        REAL(wp)     ::  xi       !<
        REAL(wp)     ::  eta      !<
@@ -3122 +3468 @@
        ENDIF
 !
-!--    First determine kcto and kctw that are the coarse-grid upper bounds for
-!--    index k
+!--    First determine kcto and kctw which refer to the uppermost 
+!--    coarse-grid levels below the child top-boundary level.
        kk = 0
        DO  WHILE ( cg%zu(kk) <= zu(nzt) )
@@ -3135 +3481 @@
        ENDDO
        kctw = kk - 1
-
-       ALLOCATE( iflu(icl:icr) )
-       ALLOCATE( iflo(icl:icr) )
-       ALLOCATE( ifuu(icl:icr) )
-       ALLOCATE( ifuo(icl:icr) )
-       ALLOCATE( jflv(jcs:jcn) )
-       ALLOCATE( jflo(jcs:jcn) )
-       ALLOCATE( jfuv(jcs:jcn) )
-       ALLOCATE( jfuo(jcs:jcn) )
+!AH
+       write(9,"('kcto, kctw = ', 2(i3,2x))") kcto, kctw
+
+!AH
+!       ALLOCATE( iflu(icl:icr) )
+!       ALLOCATE( iflo(icl:icr) )
+!       ALLOCATE( ifuu(icl:icr) )
+!       ALLOCATE( ifuo(icl:icr) )
+!       ALLOCATE( jflv(jcs:jcn) )
+!       ALLOCATE( jflo(jcs:jcn) )
+!       ALLOCATE( jfuv(jcs:jcn) )
+!       ALLOCATE( jfuo(jcs:jcn) )
+!       
+!       ALLOCATE( iflu(icl-1:icr+1) )
+!       ALLOCATE( iflo(icl-1:icr+1) )
+!       ALLOCATE( ifuu(icl-1:icr+1) )
+!       ALLOCATE( ifuo(icl-1:icr+1) )
+!       ALLOCATE( jflv(jcs-1:jcn+1) )
+!       ALLOCATE( jflo(jcs-1:jcn+1) )
+!       ALLOCATE( jfuv(jcs-1:jcn+1) )
+!       ALLOCATE( jfuo(jcs-1:jcn+1) )
+!
+       icla = coarse_bound_aux(1)
+       icra = coarse_bound_aux(2)
+       jcsa = coarse_bound_aux(3)
+       jcna = coarse_bound_aux(4)
+       ALLOCATE( iflu(icla:icra) )
+       ALLOCATE( iflo(icla:icra) )
+       ALLOCATE( ifuu(icla:icra) )
+       ALLOCATE( ifuo(icla:icra) )
+       ALLOCATE( jflv(jcsa:jcna) )
+       ALLOCATE( jflo(jcsa:jcna) )
+       ALLOCATE( jfuv(jcsa:jcna) )
+       ALLOCATE( jfuo(jcsa:jcna) )       
+!AH
        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) )
+!AH
+!       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) )
+!
+!       ALLOCATE( ijkfc_u(0:cg%nz+1,jcs-1:jcn+1,icl-1:icr+1) )
+!       ALLOCATE( ijkfc_v(0:cg%nz+1,jcs-1:jcn+1,icl-1:icr+1) )
+!       ALLOCATE( ijkfc_w(0:cg%nz+1,jcs-1:jcn+1,icl-1:icr+1) )
+!       ALLOCATE( ijkfc_s(0:cg%nz+1,jcs-1:jcn+1,icl-1:icr+1) )
+!
+       ALLOCATE( ijkfc_u(0:cg%nz+1,jcsa:jcna,icla:icra) )
+       ALLOCATE( ijkfc_v(0:cg%nz+1,jcsa:jcna,icla:icra) )
+       ALLOCATE( ijkfc_w(0:cg%nz+1,jcsa:jcna,icla:icra) )
+       ALLOCATE( ijkfc_s(0:cg%nz+1,jcsa:jcna,icla:icra) )
+!AH
 
        ijkfc_u = 0
@@ -3163 +3546 @@
        tolerance = 0.000001_wp * dx
        istart = nxlg
-       DO  ii = icl, icr-1
+!AH       DO  ii = icl, icr
+!AH       DO  ii = icl-1, icr+1
+          DO  ii = icla, icra
+
 !
 !--       In case the child and parent grid lines match in x 
@@ -3196 +3582 @@
              istart = iflu(ii)
           ENDIF
-!AH
+!
+!--       Print out the index bounds for checking and debugging purposes
           write(9,"('pmci_init_anterp_tophat, ii, iflu, ifuu: ', 3(i4,2x))")    &
                ii, iflu(ii), ifuu(ii)
@@ -3202 +3589 @@
 
        ENDDO
-       iflu(icr) = nxrg
-       ifuu(icr) = nxrg
+       write(9,*)
 !
 !--    i-indices of others for each ii-index value
 !--    ii=icr is redundant for anterpolation
        istart = nxlg
-       DO  ii = icl, icr-1
+!AH       DO  ii = icl, icr
+!AH       DO  ii = icl-1, icr+1
+       DO  ii = icla, icra
           i = istart
           DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) )  .AND.     &
@@ -3223 +3611 @@
           ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg )
           istart = iflo(ii)
+!
+!--       Print out the index bounds for checking and debugging purposes
+          write(9,"('pmci_init_anterp_tophat, ii, iflo, ifuo: ', 3(i4,2x))")    &
+               ii, iflo(ii), ifuo(ii)
+          flush(9)
        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
+       write(9,*)
 !
 !--    j-indices of v for each jj-index value
@@ -3236 +3623 @@
        tolerance = 0.000001_wp * dy
        jstart = nysg
-       DO  jj = jcs, jcn-1
+!AH       DO  jj = jcs, jcn
+!AH       DO  jj = jcs-1, jcn+1
+       DO  jj = jcsa, jcna
 !
 !--       In case the child and parent grid lines match in y 
@@ -3269 +3658 @@
              jstart = jflv(jj)
           ENDIF
-!AH
+!
+!--       Print out the index bounds for checking and debugging purposes
           write(9,"('pmci_init_anterp_tophat, jj, jflv, jfuv: ', 3(i4,2x))")    &
                jj, jflv(jj), jfuv(jj)
           flush(9)
-
        ENDDO
-       jflv(jcn) = nyng
-       jfuv(jcn) = nyng
+       write(9,*)
 !
 !--    j-indices of others for each jj-index value
 !--    jj=jcn is redundant for anterpolation
        jstart = nysg
-       DO  jj = jcs, jcn-1
+!AH       DO  jj = jcs, jcn
+!AH       DO  jj = jcs-1, jcn+1
+       DO  jj = jcsa, jcna
           j = jstart
           DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) )  .AND.     &
@@ -3296 +3686 @@
           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))")    &
+!
+!--       Print out the index bounds for checking and debugging purposes
+          write(9,"('pmci_init_anterp_tophat, jj, jflo, jfuo: ', 3(i4,2x))")    &
                jj, jflo(jj), jfuo(jj)
           flush(9)
-
        ENDDO
-       jflo(jcn) = nyng
-       jfuo(jcn) = nyng
+       write(9,*)
 !
 !--    k-indices of w for each kk-index value
@@ -3347 +3736 @@
           ENDIF
 !AH
-          write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 3(i4,2x))") &
-               kk, kflw(kk), kfuw(kk)
+          write(9,"('pmci_init_anterp_tophat, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))") &
+               kk, kflw(kk), kfuw(kk), nzt,  cg%zu(kk), cg%zw(kk)
           flush(9)
-
        ENDDO
+       write(9,*)
 !
 !--    k-indices of others for each kk-index value
@@ -3360 +3749 @@
 !--    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
-!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 )
-!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
@@ -3377 +3761 @@
           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)
+       ENDDO
+!
+!--    Set the k-index bounds separately for the parent-grid cells cg%nz and cg%nz+1. 
+!--    Index bounds for cg%nz are needed for the reversibility correction.
+       kflo(cg%nz)   = nzt+1    ! Needed for the reversibility correction
+       kfuo(cg%nz)   = nzt+kgsr ! Needed for the reversibility correction
+       kflo(cg%nz+1) = nzt+kgsr ! Obsolete
+       kfuo(cg%nz+1) = nzt+kgsr ! Obsolete 
+
+       DO  kk = 1, cg%nz+1
+          write(9,"('pmci_init_anterp_tophat, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))") &
+               kk, kflo(kk), kfuo(kk), nzt,  cg%zu(kk), cg%zw(kk)
           flush(9)
        ENDDO
+       write(9,*)
+!AH
+
 !
 !--    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
-!AH             DO kk = 0, kcto+1
+!--    This information is needed in anterpolation and in reversibility 
+!--    correction in interpolation. For the reversibility correction, ijkfc-
+!--    information is needed also beyond the indices for wall_flags_0-masking 
+!--    must be modified here in the boundary-normal direction (iw, jw, kw). 
+!AH       DO  ii = icl, icr
+!AH          DO  jj = jcs, jcn
+!AH       DO  ii = icl-1, icr+1
+!AH          DO  jj = jcs-1, jcn+1
+       DO  ii = icla, icra
+          DO  jj = jcsa, jcna
              DO kk = 0, cg%nz+1
 !
 !--             u-component
                 DO  i = iflu(ii), ifuu(ii)
+                   iw = MAX( MIN( i, nx+1 ), -1 )
                    DO  j = jflo(jj), jfuo(jj)
+                      jw = MAX( MIN( j, ny+1 ), -1 )
                       DO  k = kflo(kk), kfuo(kk)
-                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii) + MERGE( 1, 0,  &
-                                            BTEST( wall_flags_0(k,j,i), 1 ) )
+                         kw = MIN( k, nzt+1 )                
+                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii)                  &
+                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 1 ) )
                       ENDDO
                    ENDDO
@@ -3402 +3808 @@
 !--             v-component 
                 DO  i = iflo(ii), ifuo(ii)
+                   iw = MAX( MIN( i, nx+1 ), -1 )
                    DO  j = jflv(jj), jfuv(jj)
+                      jw = MAX( MIN( j, ny+1 ), -1 )
                       DO  k = kflo(kk), kfuo(kk)
-                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii) + MERGE( 1, 0,  &
-                                            BTEST( wall_flags_0(k,j,i), 2 ) )
+                         kw = MIN( k, nzt+1 )                                        
+                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii)                  &
+                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 2 ) )
                       ENDDO
                    ENDDO
@@ -3412 +3821 @@
 !--             scalars
                 DO  i = iflo(ii), ifuo(ii)
+                   iw = MAX( MIN( i, nx+1 ), -1 )
                    DO  j = jflo(jj), jfuo(jj)
+                      jw = MAX( MIN( j, ny+1 ), -1 )
                       DO  k = kflo(kk), kfuo(kk)
-                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii) + MERGE( 1, 0,  &
-                                            BTEST( wall_flags_0(k,j,i), 0 ) )
+                         kw = MIN( k, nzt+1 )
+                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii)                  &
+                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 0 ) )
                       ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
-
-!AH             DO kk = 0, kctw+1
-             DO kk = 0, cg%nz+1
 !
 !--             w-component 
                 DO  i = iflo(ii), ifuo(ii)
+                   iw = MAX( MIN( i, nx+1 ), -1 )
                    DO  j = jflo(jj), jfuo(jj)
+                      jw = MAX( MIN( j, ny+1 ), -1 )
                       DO  k = kflw(kk), kfuw(kk)
-                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii) + MERGE( 1, 0,  &
-                                            BTEST( wall_flags_0(k,j,i), 3 ) )
+                         kw = MIN( k, nzt+1 )
+                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii) + MERGE( 1, 0,   &
+                              BTEST( wall_flags_0(kw,jw,iw), 3 ) )
                       ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
-        
-          ENDDO
-       ENDDO
+
+             ENDDO  ! kk        
+          ENDDO  ! jj
+       ENDDO  ! ii
 !
 !--    Spatial under-relaxation coefficients
@@ -3445 +3856 @@
        fray(jcs:jcn) = 1.0_wp
 
-       IF ( nesting_mode /= 'vertical' )  THEN
-          DO  ii = icl, icr
-             IF ( ifuu(ii) < ( nx + 1 ) / 2 )  THEN   
-                xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) -                         &
-                     lower_left_coord_x ) ) / anterp_relax_length_l )**4
-                frax(ii) = xi / ( 1.0_wp + xi )
-             ELSE
-                xi = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx -   &
-                                      cg%coord_x(ii) ) ) /                     &
-                       anterp_relax_length_r )**4
-                frax(ii) = xi / ( 1.0_wp + xi )                
-             ENDIF
-          ENDDO
-
-          DO  jj = jcs, jcn
-             IF ( jfuv(jj) < ( ny + 1 ) / 2 )  THEN 
-                eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) -                        &
-                     lower_left_coord_y ) ) / anterp_relax_length_s )**4
-                fray(jj) = eta / ( 1.0_wp + eta )
-             ELSE
-                eta = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy -  &
-                                       cg%coord_y(jj)) ) /                     &
-                        anterp_relax_length_n )**4
-                fray(jj) = eta / ( 1.0_wp + eta )
-             ENDIF
-          ENDDO
-       ENDIF
+!AH       IF ( nesting_mode /= 'vertical' )  THEN
+!AH          DO  ii = icl, icr
+!AH             IF ( ifuu(ii) < ( nx + 1 ) / 2 )  THEN   
+!AH                xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) -                         &
+!AH                     lower_left_coord_x ) ) / anterp_relax_length_l )**4
+!AH                frax(ii) = xi / ( 1.0_wp + xi )
+!AH             ELSE
+!AH                xi = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx -   &
+!AH                                      cg%coord_x(ii) ) ) /                     &
+!AH                       anterp_relax_length_r )**4
+!AH                frax(ii) = xi / ( 1.0_wp + xi )                
+!AH             ENDIF
+!AH          ENDDO
+!AH
+!AH          DO  jj = jcs, jcn
+!AH             IF ( jfuv(jj) < ( ny + 1 ) / 2 )  THEN 
+!AH                eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) -                        &
+!AH                     lower_left_coord_y ) ) / anterp_relax_length_s )**4
+!AH                fray(jj) = eta / ( 1.0_wp + eta )
+!AH             ELSE
+!AH                eta = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy -  &
+!AH                                       cg%coord_y(jj)) ) /                     &
+!AH                        anterp_relax_length_n )**4
+!AH                fray(jj) = eta / ( 1.0_wp + eta )
+!AH             ENDIF
+!AH          ENDDO
+!AH       ENDIF
       
        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 )
-       ENDDO
+       fraz(0:kcto) = 1.0_wp
+!AH       DO  kk = 0, kcto
+!AH          zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4
+!AH          fraz(kk) = zeta / ( 1.0_wp + zeta )
+!AH       ENDDO
 
     END SUBROUTINE pmci_init_anterp_tophat
-
-
-
-    SUBROUTINE pmci_init_tkefactor
-
-!
-!--    Computes the scaling factor for the SGS TKE from coarse grid to be used
-!--    as BC for the fine grid. Based on the Kolmogorov energy spectrum
-!--    for the inertial subrange and assumption of sharp cut-off of the resolved
-!--    energy spectrum. Near the surface, the reduction of TKE is made
-!--    smaller than further away from the surface.
-!--    Please note, in case parent and child model operate in RANS mode, 
-!--    TKE is not grid depenedent and weighting factor is one. 
-
-       IMPLICIT NONE
-
-       INTEGER(iwp)        ::  k                     !< index variable along z
-       INTEGER(iwp)        ::  k_wall                !< topography-top index along z
-       INTEGER(iwp)        ::  kc                    !<
-
-       REAL(wp), PARAMETER ::  cfw = 0.2_wp          !<
-       REAL(wp), PARAMETER ::  c_tkef = 0.6_wp       !<
-       REAL(wp)            ::  fw                    !<
-       REAL(wp), PARAMETER ::  fw0 = 0.9_wp          !<
-       REAL(wp)            ::  glsf                  !<
-       REAL(wp)            ::  glsc                  !<
-       REAL(wp)            ::  height                !<
-       REAL(wp), PARAMETER ::  p13 = 1.0_wp/3.0_wp   !<
-       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !<        
-
-!
-       IF ( .NOT. rans_mode  .AND.  .NOT. rans_mode_parent )  THEN 
-          IF ( bc_dirichlet_l )  THEN
-             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
-             tkefactor_l = 0.0_wp
-             i = nxl - 1
-             DO  j = nysg, nyng
-                k_wall = get_topography_top_index_ji( j, i, 's' )
-
-                DO  k = k_wall + 1, nzt
-                   kc     = kco(k) + 1
-                   glsf   = ( dx * dy * dzu(k) )**p13
-                   glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
-                   height = zu(k) - zu(k_wall)
-                   fw     = EXP( -cfw * height / glsf )
-                   tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
-                                                 ( glsf / glsc )**p23 )
-                ENDDO
-                tkefactor_l(k_wall,j) = c_tkef * fw0
-             ENDDO
-          ENDIF
-
-          IF ( bc_dirichlet_r )  THEN
-             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
-             tkefactor_r = 0.0_wp
-             i = nxr + 1
-             DO  j = nysg, nyng
-                k_wall = get_topography_top_index_ji( j, i, 's' )
-
-                DO  k = k_wall + 1, nzt
-                   kc     = kco(k) + 1
-                   glsf   = ( dx * dy * dzu(k) )**p13
-                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                   height = zu(k) - zu(k_wall)
-                   fw     = EXP( -cfw * height / glsf )
-                   tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
-                                                 ( glsf / glsc )**p23 )
-                ENDDO
-                tkefactor_r(k_wall,j) = c_tkef * fw0
-             ENDDO
-          ENDIF
-
-          IF ( bc_dirichlet_s )  THEN
-             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
-             tkefactor_s = 0.0_wp
-             j = nys - 1
-             DO  i = nxlg, nxrg
-                k_wall = get_topography_top_index_ji( j, i, 's' )
-                
-                DO  k = k_wall + 1, nzt
-    
-                   kc     = kco(k) + 1
-                   glsf   = ( dx * dy * dzu(k) )**p13
-                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
-                   height = zu(k) - zu(k_wall)
-                   fw     = EXP( -cfw*height / glsf )
-                   tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
-                        ( glsf / glsc )**p23 )
-                ENDDO
-                tkefactor_s(k_wall,i) = c_tkef * fw0
-             ENDDO
-          ENDIF
-
-          IF ( bc_dirichlet_n )  THEN
-             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
-             tkefactor_n = 0.0_wp
-             j = nyn + 1
-             DO  i = nxlg, nxrg
-                k_wall = get_topography_top_index_ji( j, i, 's' )
-
-                DO  k = k_wall + 1, nzt
-
-                   kc     = kco(k) + 1
-                   glsf   = ( dx * dy * dzu(k) )**p13
-                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                   height = zu(k) - zu(k_wall)
-                   fw     = EXP( -cfw * height / glsf )
-                   tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
-                                                 ( glsf / glsc )**p23 )
-                ENDDO
-                tkefactor_n(k_wall,i) = c_tkef * fw0
-             ENDDO
-          ENDIF
-
-          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
-          k = nzt
-
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-!
-!--             Determine vertical index for local topography top
-                k_wall = get_topography_top_index_ji( j, i, 's' )
-
-                kc     = kco(k) + 1
-                glsf   = ( dx * dy * dzu(k) )**p13
-                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-                height = zu(k) - zu(k_wall)
-                fw     = EXP( -cfw * height / glsf )
-                tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
-                                              ( glsf / glsc )**p23 )
-             ENDDO
-          ENDDO
-!
-!--    RANS mode
-       ELSE
-          IF ( bc_dirichlet_l )  THEN
-             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
-             tkefactor_l = 1.0_wp
-          ENDIF
-          IF ( bc_dirichlet_r )  THEN
-             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
-             tkefactor_r = 1.0_wp
-          ENDIF
-          IF ( bc_dirichlet_s )  THEN
-             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
-             tkefactor_s = 1.0_wp
-          ENDIF
-          IF ( bc_dirichlet_n )  THEN
-             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
-             tkefactor_n = 1.0_wp
-          ENDIF
-
-          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
-          tkefactor_t = 1.0_wp
-
-       ENDIF
-      
-    END SUBROUTINE pmci_init_tkefactor
 
 #endif
@@ -3667 +3921 @@
 #endif
  END SUBROUTINE pmci_setup_coordinates
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -4029 +4283 @@
     INTEGER(iwp) ::  i          !<
     INTEGER(iwp) ::  icl        !<
+    INTEGER(iwp) ::  icla       !<
+    INTEGER(iwp) ::  iclw       !<
     INTEGER(iwp) ::  icr        !<
+    INTEGER(iwp) ::  icra       !<
+    INTEGER(iwp) ::  icrw       !<
     INTEGER(iwp) ::  j          !<
     INTEGER(iwp) ::  jcn        !<
+    INTEGER(iwp) ::  jcna       !<
+    INTEGER(iwp) ::  jcnw       !<
     INTEGER(iwp) ::  jcs        !<
+    INTEGER(iwp) ::  jcsa       !<
+    INTEGER(iwp) ::  jcsw       !<
     INTEGER(iwp) ::  k          !< 
     INTEGER(iwp) ::  n          !< running index for chemical species
@@ -4043 +4305 @@
 !
 !--    Child domain boundaries in the parent index space
-       icl = coarse_bound(1)
-       icr = coarse_bound(2)
-       jcs = coarse_bound(3)
-       jcn = coarse_bound(4)
+       icl  = coarse_bound(1)
+       icr  = coarse_bound(2)
+       jcs  = coarse_bound(3)
+       jcn  = coarse_bound(4)
+       icla = coarse_bound_aux(1)
+       icra = coarse_bound_aux(2)
+       jcsa = coarse_bound_aux(3)
+       jcna = coarse_bound_aux(4)
+       iclw = coarse_bound_w(1)
+       icrw = coarse_bound_w(2)
+       jcsw = coarse_bound_w(3)
+       jcnw = coarse_bound_w(4)
+
 !
 !--    Get data from the parent
@@ -4052 +4323 @@
 !
 !--    The interpolation.
-       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
+       CALL pmci_interp_1sto_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
                                    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,   &
+       CALL pmci_interp_1sto_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,   &
                                    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,   &
+       CALL pmci_interp_1sto_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,   &
                                    r2yo, r1zw, r2zw, kctw, iflo, ifuo,         &
                                    jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' )
@@ -4065 +4336 @@
             (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.           &
                .NOT. constant_diffusion ) )  THEN
-          CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
+          CALL pmci_interp_1sto_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
                                       r2yo, r1zo, r2zo, kcto, iflo, ifuo,       &
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' )
@@ -4071 +4342 @@
 
        IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-          CALL pmci_interp_tril_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
+          CALL pmci_interp_1sto_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
                                       r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,&
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -4077 +4348 @@
 
        IF ( .NOT. neutral )  THEN
-          CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
+          CALL pmci_interp_1sto_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
                                       r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,&
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -4084 +4355 @@
        IF ( humidity )  THEN
 
-          CALL pmci_interp_tril_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, &
+          CALL pmci_interp_1sto_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, &
                                       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,   &
+             CALL pmci_interp_1sto_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,   &
                                          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,   &
+             CALL pmci_interp_1sto_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,   &
                                          r1yo, r2yo, r1zo, r2zo, kcto,         &
                                          iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
@@ -4100 +4371 @@
 
           IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-             CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,   &
+             CALL pmci_interp_1sto_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,   &
                                          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,   &
+             CALL pmci_interp_1sto_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,   &
                                          r1yo, r2yo, r1zo, r2zo, kcto,         &
                                          iflo, ifuo, jflo, jfuo, kflo, kfuo,   &
@@ -4113 +4384 @@
 
        IF ( passive_scalar )  THEN
-          CALL pmci_interp_tril_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+          CALL pmci_interp_1sto_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,  &
                                       r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -4120 +4391 @@
        IF ( air_chemistry  .AND.  nest_chemistry )  THEN
           DO  n = 1, nspec
-             CALL pmci_interp_tril_all ( chem_species(n)%conc,                 &
+             CALL pmci_interp_1sto_all ( chem_species(n)%conc,                 &
                                          chem_spec_c(:,:,:,n),                 &
                                          ico, jco, kco, r1xo, r2xo, r1yo,      &
@@ -4162 +4433 @@
 
 
-    SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
+    SUBROUTINE pmci_interp_1sto_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
                                      r1z, r2z, kct, ifl, ifu, jfl, jfu,        &
                                      kfl, kfu, ijkfc, var )
@@ -4173 +4444 @@
        CHARACTER(LEN=1), INTENT(IN) :: var  !<
 
-       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !<
-       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !<
-       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
+       INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)         ::  ic    !<
+       INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)         ::  jc    !<
+!AH       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
+       INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), 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   !<
+       REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  :: r1x   !<
+       REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  :: r2x   !<
+       REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  :: r1y   !<
+       REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  :: r2y   !<
+       REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) :: r1z   !<
+       REAL(wp), DIMENSION(nzb:nzt+kgsr), 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
+!AH       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+!AH       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+!AH       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(jcsa:jcna), 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
+!AH
+!       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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
 
        INTEGER(iwp) ::  i        !<
-       INTEGER(iwp) ::  ib       !<
-       INTEGER(iwp) ::  ie       !<
-       INTEGER(iwp) ::  ijk      !<
        INTEGER(iwp) ::  j        !<
-       INTEGER(iwp) ::  jb       !<
-       INTEGER(iwp) ::  je       !<
        INTEGER(iwp) ::  k        !<
-       INTEGER(iwp) ::  k_wall   !<
-       INTEGER(iwp) ::  k1       !<
-       INTEGER(iwp) ::  kb       !<
-       INTEGER(iwp) ::  kbc      !<
        INTEGER(iwp) ::  l        !<
+       INTEGER(iwp) ::  lb       !<
+       INTEGER(iwp) ::  le       !<
        INTEGER(iwp) ::  m        !<
+       INTEGER(iwp) ::  mb       !<
+       INTEGER(iwp) ::  me       !<
        INTEGER(iwp) ::  n        !<
        INTEGER(iwp) ::  var_flag !<
 
-       REAL(wp) ::  cellsum    !<
-       REAL(wp) ::  cellsumd   !<
-       REAL(wp) ::  fk         !<
-       REAL(wp) ::  fkj        !<
-       REAL(wp) ::  fkjp       !<
-       REAL(wp) ::  fkp        !<
-       REAL(wp) ::  fkpj       !<
-       REAL(wp) ::  fkpjp      !<
-       REAL(wp) ::  logratio   !<
-       REAL(wp) ::  logzuc1    !<
-       REAL(wp) ::  rcorr      !<
-       REAL(wp) ::  rcorr_ijk  !<       
-       REAL(wp) ::  zuc1       !<
-       REAL(wp) ::  z0_topo    !<  roughness at vertical walls
-
-
-       ib = nxl
-       ie = nxr
-       jb = nys
-       je = nyn
-       kb = 0
+
+       lb = icl
+       le = icr
+       mb = jcs
+       me = jcn
+
        IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
-             ib = nxl - 1
+             lb = icl + 1
 !
 !--          For u, nxl is a ghost node, but not for the other variables
              IF ( var == 'u' )  THEN
-                ib = nxl
+                lb = icl + 2
              ENDIF
           ENDIF
           IF ( bc_dirichlet_s )  THEN
-             jb = nys - 1
+             mb = jcs + 1
 !
 !--          For v, nys is a ghost node, but not for the other variables
              IF ( var == 'v' )  THEN
-                jb = nys
+                mb = jcs + 2
              ENDIF
           ENDIF
           IF ( bc_dirichlet_r )  THEN
-             ie = nxr + 1
+             le = icr - 1
           ENDIF
           IF ( bc_dirichlet_n )  THEN
-             je = nyn + 1
+             me = jcn - 1
           ENDIF
        ENDIF
-
+!
+!--    Is masking needed here, think about it.
        IF ( var == 'u' )  THEN 
           var_flag = 1 
@@ -4266 +4526 @@
        ELSE
           var_flag = 0
-       ENDIF
-!
-!--    Trilinear interpolation.
-       DO  i = ib, ie
-          DO  j = jb, je
-!
-!--          Determine the vertical index of the first node above the
-!--          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             
-             DO  k = kb, nzt + 1
-                l = ic(i)
-                m = jc(j)
-                n = kc(k)
-                fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
-                fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
-                fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
-                fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
-                fk       = r1y(j) * fkj  + r2y(j) * fkjp
-                fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
-                f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+       ENDIF       
+
+       f(:,:,:) = 0.0_wp
+
+       IF  ( var == 'u' )  THEN
+
+          DO  l = lb, le
+             DO  m = mb, me
+                DO n = 0, kct 
+
+                   DO  i = ifl(l), ifl(l+1)-1
+                      DO  j = jfl(m), jfu(m)
+                         DO  k = kfl(n), kfu(n)
+                            f(k,j,i) = fc(n,m,l)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   
+                ENDDO
              ENDDO
           ENDDO
-       ENDDO
-!
-!--    Correct the interpolated values of u and v in near-wall nodes, i.e. in
-!--    the nodes below the coarse-grid nodes with k=1. The corrction is only
-!--    made over horizontal wall surfaces in this phase. For the nest boundary
-!--    conditions, a corresponding correction is made for all vertical walls,
-!--    too.
-       IF ( constant_flux_layer .AND. ( var == 'u' .OR. var == 'v' ) )  THEN
-          z0_topo = roughness_length
-          DO  i = ib, nxr
-             DO  j = jb, nyn
-!
-!--             Determine vertical index of topography top at grid point (j,i)
-                k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
-!
-!--             kbc is the first coarse-grid point above the surface
-                kbc = 1
-                DO  WHILE ( cg%zu(kbc) < zu(k_wall) )
-                   kbc = kbc + 1
+
+       ELSE IF  ( var == 'v' )  THEN
+
+          DO  l = lb, le
+             DO  m = mb, me
+                DO n = 0, kct 
+               
+                   DO i = ifl(l), ifu(l)
+                      DO  j = jfl(m), jfl(m+1)-1
+                         DO  k = kfl(n), kfu(n)
+                            f(k,j,i) = fc(n,m,l) 
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
                 ENDDO
-                zuc1 = cg%zu(kbc)
-                k1   = k_wall + 1
-                DO  WHILE ( zu(k1) < zuc1 )
-                   k1 = k1 + 1
-                ENDDO
-                logzuc1 = LOG( ( zu(k1) - zu(k_wall) ) / z0_topo )
-
-                k = k_wall + 1
-                DO  WHILE ( zu(k) < zuc1 )
-                   logratio = ( LOG( ( zu(k) - zu(k_wall) ) / z0_topo ) ) /    &
-                                logzuc1
-                   f(k,j,i) = logratio * f(k1,j,i)
-                   k  = k + 1
-                ENDDO
-                f(k_wall,j,i) = 0.0_wp
              ENDDO
           ENDDO
 
-       ELSEIF ( var == 'w' )  THEN
-
-          DO  i = ib, nxr
-              DO  j = jb, nyn
-!
-!--              Determine vertical index of topography top at grid point (j,i)
-                 k_wall = get_topography_top_index_ji( j, i, 'w' )
-
-                 f(k_wall,j,i) = 0.0_wp
-              ENDDO
-           ENDDO
-
-       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
+       ELSE IF  ( var == 'w' )  THEN
+
+          DO  l = lb, le
+             DO  m = mb, me
+                DO n = 1, kct + 1   ! It is important to go up to kct+1  
+
+                   DO i = ifl(l), ifu(l)
+                      DO  j = jfl(m), jfu(m)
+                         f(nzb,j,i) = 0.0_wp   ! Because the n-loop starts from n=1 instead of 0
+                         DO  k = kfu(n-1)+1, kfu(n) 
+                            f(k,j,i) = fc(n,m,l)
+                         ENDDO
                       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
+
+       ELSE   ! scalars
+
+          DO  l = lb, le
+             DO  m = mb, me
+                DO n = 0, kct 
+                   
+                   DO i = ifl(l), ifu(l)
+                      DO  j = jfl(m), jfu(m)
+                         DO  k = kfl(n), kfu(n)
+                            f(k,j,i) = fc(n,m,l)
+                         ENDDO
                       ENDDO
                    ENDDO
+                   
                 ENDDO
-             
-             ENDDO  ! n
-          ENDDO  ! m
-       ENDDO  ! l
-
-    END SUBROUTINE pmci_interp_tril_all
+             ENDDO
+          ENDDO
+
+       ENDIF  ! var
+
+    END SUBROUTINE pmci_interp_1sto_all
 
 #endif
@@ -4674 +4878 @@
     IMPLICIT NONE
 
-    INTEGER(iwp), INTENT(IN) ::  direction   !<
+    INTEGER(iwp), INTENT(IN) ::  direction  !< Transfer direction: parent_to_child or child_to_parent
 
 #if defined( __parallel )
-    INTEGER(iwp) ::  icl         !<
-    INTEGER(iwp) ::  icr         !<
-    INTEGER(iwp) ::  jcs         !<
-    INTEGER(iwp) ::  jcn         !<
-    
-    REAL(wp), DIMENSION(1) ::  dtl         !<
+    INTEGER(iwp) ::  icl         !< Parent-grid array index bound, left
+    INTEGER(iwp) ::  icr         !< Parent-grid array index bound, right
+    INTEGER(iwp) ::  jcn         !< Parent-grid array index bound, north
+    INTEGER(iwp) ::  jcs         !< Parent-grid array index bound, south
+    INTEGER(iwp) ::  icla        !< Auxiliary-array (index-mapping etc) index bound, left
+    INTEGER(iwp) ::  icra        !< Auxiliary-array (index-mapping etc) index bound, right
+    INTEGER(iwp) ::  jcna        !< Auxiliary-array (index-mapping etc) index bound, north
+    INTEGER(iwp) ::  jcsa        !< Auxiliary-array (index-mapping etc) index bound, south
+    INTEGER(iwp) ::  iclw        !< Parent-grid work array index bound, left
+    INTEGER(iwp) ::  icrw        !< Parent-grid work array index bound, right
+    INTEGER(iwp) ::  jcnw        !< Parent-grid work array index bound, north
+    INTEGER(iwp) ::  jcsw        !< Parent-grid work array index bound, south
+
+    REAL(wp), DIMENSION(1) ::  dtl         !< Time step size
 
 
@@ -4689 +4901 @@
 !
 !--    Child domain boundaries in the parent indice space.
-       icl = coarse_bound(1)
-       icr = coarse_bound(2)
-       jcs = coarse_bound(3)
-       jcn = coarse_bound(4)
+       icl  = coarse_bound(1)
+       icr  = coarse_bound(2)
+       jcs  = coarse_bound(3)
+       jcn  = coarse_bound(4)
+       icla = coarse_bound_aux(1)
+       icra = coarse_bound_aux(2)
+       jcsa = coarse_bound_aux(3)
+       jcna = coarse_bound_aux(4)
+       iclw = coarse_bound_w(1)
+       icrw = coarse_bound_w(2)
+       jcsw = coarse_bound_w(3)
+       jcnw = coarse_bound_w(4)
 
        IF ( direction == parent_to_child )  THEN
@@ -4736 +4956 @@
 !--    horizontal boundaries
        IF ( nesting_mode /= 'vertical' )  THEN
-       
 !
 !--       Left border pe:
           IF ( bc_dirichlet_l )  THEN
 
-             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
+             CALL pmci_interp_1sto_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,      &
                                        r1yo, r2yo, r1zo, r2zo,                  &
                                        logc_u_l, logc_ratio_u_l,                &
@@ -4748 +4967 @@
                                        kfuo, ijkfc_u, 'l', 'u' )
 
-             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
+             CALL pmci_interp_1sto_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,      &
                                        r1yv, r2yv, r1zo, r2zo,                  &
                                        logc_v_l, logc_ratio_v_l,                &
@@ -4755 +4974 @@
                                        kfuo, ijkfc_v, 'l', 'v' )
 
-             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
+             CALL pmci_interp_1sto_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,      &
                                        r1yo, r2yo, r1zw, r2zw,                  &
                                        logc_w_l, logc_ratio_w_l,                &
@@ -4765 +4984 @@
                   (  .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,  &
+!                CALL pmci_interp_1sto_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
 !                                          r1yo, r2yo, r1zo, r2zo,              &
 !                                          logc_w_l, logc_ratio_w_l,            &
@@ -4780 +4999 @@
 
              IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
+                CALL pmci_interp_1sto_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_l, logc_ratio_w_l,            &
@@ -4789 +5008 @@
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
+                CALL pmci_interp_1sto_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_l, logc_ratio_w_l,            &
@@ -4799 +5018 @@
              IF ( humidity )  THEN
 
-                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
+                CALL pmci_interp_1sto_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_l, logc_ratio_w_l,            &
@@ -4807 +5026 @@
 
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
-                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( qc, qcc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_l, logc_ratio_w_l,         &
@@ -4815 +5034 @@
                                              kflo, kfuo, ijkfc_s, 'l', 's' )  
 
-                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( nc, ncc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_l, logc_ratio_w_l,         &
@@ -4825 +5044 @@
 
                 IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
-                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_l, logc_ratio_w_l,         &
@@ -4833 +5052 @@
                                              kflo, kfuo, ijkfc_s, 'l', 's' ) 
 
-                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_l, logc_ratio_w_l,         &
@@ -4845 +5064 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
+                CALL pmci_interp_1sto_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_l, logc_ratio_w_l,            &
@@ -4856 +5075 @@
              IF ( air_chemistry  .AND.  nest_chemistry )  THEN
                 DO  n = 1, nspec
-                   CALL pmci_interp_tril_lr( chem_species(n)%conc,             &
+                   CALL pmci_interp_1sto_lr( chem_species(n)%conc,             &
                                              chem_spec_c(:,:,:,n),             &
                                              ico, jco, kco, r1xo, r2xo,        &
@@ -4873 +5092 @@
           IF ( bc_dirichlet_r )  THEN
              
-             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
+             CALL pmci_interp_1sto_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_r, logc_ratio_u_r,               &
@@ -4880 +5099 @@
                                        kfuo, ijkfc_u, 'r', 'u' )
 
-             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
+             CALL pmci_interp_1sto_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
                                        r1yv, r2yv, r1zo, r2zo,                 &
                                        logc_v_r, logc_ratio_v_r,               &
@@ -4887 +5106 @@
                                        kfuo, ijkfc_v, 'r', 'v' )
 
-             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
+             CALL pmci_interp_1sto_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
                                        r1yo, r2yo, r1zw, r2zw,                 &
                                        logc_w_r, logc_ratio_w_r,               &
@@ -4897 +5116 @@
                   (  .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,  &
+!                CALL pmci_interp_1sto_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
 !                                          r1yo,r2yo, r1zo, r2zo,               &
 !                                          logc_w_r, logc_ratio_w_r,            &
@@ -4911 +5130 @@
 
              IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
+                CALL pmci_interp_1sto_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
                                           r2xo, r1yo,r2yo, r1zo, r2zo,         &
                                           logc_w_r, logc_ratio_w_r,            &
@@ -4921 +5140 @@
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
+                CALL pmci_interp_1sto_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_r, logc_ratio_w_r,            &
@@ -4930 +5149 @@
 
              IF ( humidity )  THEN
-                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
+                CALL pmci_interp_1sto_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_r, logc_ratio_w_r,            &
@@ -4939 +5158 @@
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
-                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( qc, qcc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_r, logc_ratio_w_r,         &
@@ -4947 +5166 @@
                                              kflo, kfuo, ijkfc_s, 'r', 's' ) 
      
-                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( nc, ncc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_r, logc_ratio_w_r,         &
@@ -4960 +5179 @@
 
      
-                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_r, logc_ratio_w_r,         &
@@ -4969 +5188 @@
                                              'r', 's' ) 
 
-                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_lr( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_r, logc_ratio_w_r,         &
@@ -4982 +5201 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
+                CALL pmci_interp_1sto_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_r, logc_ratio_w_r,            &
@@ -4993 +5212 @@
              IF ( air_chemistry  .AND.  nest_chemistry )  THEN
                 DO  n = 1, nspec
-                   CALL pmci_interp_tril_lr( chem_species(n)%conc,             &
+                   CALL pmci_interp_1sto_lr( chem_species(n)%conc,             &
                                              chem_spec_c(:,:,:,n),             &
                                              ico, jco, kco, r1xo, r2xo,        &
@@ -5009 +5228 @@
           IF ( bc_dirichlet_s )  THEN
 
-             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
+             CALL pmci_interp_1sto_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,    &
+                                       kcto, iflo, ifuo, jflv, jfuv, kflo,     &
+                                       kfuo, ijkfc_v, 's', 'v' )
+
+             CALL pmci_interp_1sto_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,    &
+                                       kctw, iflo, ifuo, jflo, jfuo, kflw,     &
+                                       kfuw, ijkfc_w, 's','w' )
+
+             CALL pmci_interp_1sto_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_s, logc_ratio_u_s,               &
@@ -5016 +5249 @@
                                        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,    &
-                                       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,    &
-                                       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,  &
+!                CALL pmci_interp_1sto_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
 !                                          r1yo, r2yo, r1zo, r2zo,              &
 !                                          logc_w_s, logc_ratio_w_s,            &
@@ -5047 +5266 @@
 
              IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
+                CALL pmci_interp_1sto_sn( diss, dissc,  ico, jco, kco, r1xo,   &
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_s, logc_ratio_w_s,            &
@@ -5056 +5275 @@
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
+                CALL pmci_interp_1sto_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_s, logc_ratio_w_s,            &
@@ -5065 +5284 @@
 
              IF ( humidity )  THEN
-                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
+                CALL pmci_interp_1sto_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
                                           r1yo,r2yo, r1zo, r2zo,               &
                                           logc_w_s, logc_ratio_w_s,            &
@@ -5074 +5293 @@
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
-                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( qc, qcc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
                                              logc_w_s, logc_ratio_w_s,         &
@@ -5082 +5301 @@
                                              kflo, kfuo, ijkfc_s, 's', 's' )
 
-                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( nc, ncc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
                                              logc_w_s, logc_ratio_w_s,         &
@@ -5094 +5313 @@
                 IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
 
-                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
                                              logc_w_s, logc_ratio_w_s,         &
@@ -5102 +5321 @@
                                              kflo, kfuo, ijkfc_s, 's', 's' )
 
-                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
                                              logc_w_s, logc_ratio_w_s,         &
@@ -5115 +5334 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
+                CALL pmci_interp_1sto_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
                                           r1yo,r2yo, r1zo, r2zo,               &
                                           logc_w_s, logc_ratio_w_s,            &
@@ -5126 +5345 @@
              IF ( air_chemistry  .AND.  nest_chemistry )  THEN
                 DO  n = 1, nspec
-                   CALL pmci_interp_tril_sn( chem_species(n)%conc,             &
+                   CALL pmci_interp_1sto_sn( chem_species(n)%conc,             &
                                              chem_spec_c(:,:,:,n),             &
                                              ico, jco, kco, r1xo, r2xo,        &
@@ -5142 +5361 @@
           IF ( bc_dirichlet_n )  THEN
              
-             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
+             CALL pmci_interp_1sto_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
                                        r1yo, r2yo, r1zo, r2zo,                 &
                                        logc_u_n, logc_ratio_u_n,               &
@@ -5149 +5368 @@
                                        kfuo, ijkfc_u, 'n', 'u' )
 
-             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
+             CALL pmci_interp_1sto_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
                                        r1yv, r2yv, r1zo, r2zo,                 &
                                        logc_v_n, logc_ratio_v_n,               &
@@ -5156 +5375 @@
                                        kfuo, ijkfc_v, 'n', 'v' )
 
-             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
+             CALL pmci_interp_1sto_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
                                        r1yo, r2yo, r1zw, r2zw,                 &
                                        logc_w_n, logc_ratio_w_n,               &
@@ -5166 +5385 @@
                   (  .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,  &
+!                CALL pmci_interp_1sto_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
 !                                          r1yo, r2yo, r1zo, r2zo,              &
 !                                          logc_w_n, logc_ratio_w_n,            &
@@ -5180 +5399 @@
 
              IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
+                CALL pmci_interp_1sto_sn( diss, dissc,  ico, jco, kco, r1xo,   &
                                           r2xo, r1yo, r2yo, r1zo, r2zo,        &
                                           logc_w_n, logc_ratio_w_n,            &
@@ -5190 +5409 @@
 
              IF ( .NOT. neutral )  THEN
-                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
+                CALL pmci_interp_1sto_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_n, logc_ratio_w_n,            &
@@ -5199 +5418 @@
 
              IF ( humidity )  THEN
-                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
+                CALL pmci_interp_1sto_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_n, logc_ratio_w_n,            &
@@ -5208 +5427 @@
                 IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
-                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( qc, qcc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_n, logc_ratio_w_n,         &
@@ -5216 +5435 @@
                                              kflo, kfuo, ijkfc_s, 'n', 's' )
 
-                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( nc, ncc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_u_n, logc_ratio_u_n,         &
@@ -5228 +5447 @@
                 IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
 
-                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( qr, qrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_n, logc_ratio_w_n,         &
@@ -5236 +5455 @@
                                              kflo, kfuo, ijkfc_s, 'n', 's' )
 
-                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
+                   CALL pmci_interp_1sto_sn( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
                                              logc_w_n, logc_ratio_w_n,         &
@@ -5249 +5468 @@
 
              IF ( passive_scalar )  THEN
-                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
+                CALL pmci_interp_1sto_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
                                           r1yo, r2yo, r1zo, r2zo,              &
                                           logc_w_n, logc_ratio_w_n,            &
@@ -5260 +5479 @@
              IF ( air_chemistry  .AND.  nest_chemistry )  THEN
                 DO  n = 1, nspec
-                   CALL pmci_interp_tril_sn( chem_species(n)%conc,             &
+                   CALL pmci_interp_1sto_sn( chem_species(n)%conc,             &
                                              chem_spec_c(:,:,:,n),             &
                                              ico, jco, kco, r1xo, r2xo,        &
@@ -5276 +5495 @@
 !
 !--    All PEs are top-border PEs
-       CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
+       CALL pmci_interp_1sto_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
+                                r2yo, r1zw, r2zw, kctw, iflo, ifuo,            &
+                                jflo, jfuo, kflw, kfuw, ijkfc_w, 'w' )
+       CALL pmci_interp_1sto_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
                                 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,      &
+       CALL pmci_interp_1sto_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,      &
                                 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, 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,   &
+!          CALL pmci_interp_1sto_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
 !                                   r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
 !                                   jflo, jfuo, kflo, kfuo, ijkfc_s, 'e' )
@@ -5299 +5519 @@
 
        IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
-          CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
+          CALL pmci_interp_1sto_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
                                    r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,   &
                                    jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5305 +5525 @@
 
        IF ( .NOT. neutral )  THEN
-          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo,         &
+          CALL pmci_interp_1sto_t( pt, ptc, ico, jco, kco, r1xo, r2xo,         &
                                    r1yo, r2yo, r1zo, r2zo, kcto, iflo, ifuo,   &
                                    jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5312 +5532 @@
        IF ( humidity )  THEN
 
-          CALL pmci_interp_tril_t( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo,    &
+          CALL pmci_interp_1sto_t( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo,    &
                                    r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
                                    jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5318 +5538 @@
           IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
 
-             CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
+             CALL pmci_interp_1sto_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
                                       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,&
+             CALL pmci_interp_1sto_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
                                       r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5331 +5551 @@
 
 
-             CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
+             CALL pmci_interp_1sto_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
                                       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,&
+             CALL pmci_interp_1sto_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,&
                                       r2yo, r1zo, r2zo, kcto, iflo, ifuo,      &
                                       jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5344 +5564 @@
 
        IF ( passive_scalar )  THEN
-          CALL pmci_interp_tril_t( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,     &
+          CALL pmci_interp_1sto_t( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,     &
                                    r2yo, r1zo, r2zo, kcto, iflo, ifuo,         &
                                    jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
@@ -5351 +5571 @@
        IF ( air_chemistry  .AND.  nest_chemistry )  THEN
           DO  n = 1, nspec
-             CALL pmci_interp_tril_t( chem_species(n)%conc,                    &
+             CALL pmci_interp_1sto_t( chem_species(n)%conc,                    &
                                       chem_spec_c(:,:,:,n),                    &
                                       ico, jco, kco, r1xo, r2xo,               &
@@ -5444 +5664 @@
 
 
-   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
+   SUBROUTINE pmci_interp_1sto_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
                                    r2z, logc, logc_ratio, logc_kbounds,        &
                                    nzt_topo_nestbc,                            &
@@ -5464 +5684 @@
       REAL(wp), DIMENSION(1:2,0:ncorr-1,nzb:nzt_topo_nestbc,nys:nyn),          &
                                       INTENT(IN)    ::  logc_ratio   !<
-      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     !<
+!AH      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r1x     !<
+!AH      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r2x     !<
+!AH      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r1y     !<
+!AH      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r2y     !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x     !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x     !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y     !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y     !<
+
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z     !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z     !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z     !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z     !<
       
-      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic     !<
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc     !<
-      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !<
+
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)         ::  ic     !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)         ::  jc     !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), INTENT(IN)        ::  kc     !<
       INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nys:nyn),                &
                                           INTENT(IN)           :: logc   !<
@@ -5479 +5708 @@
 
       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(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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH
+
 !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
@@ -5489 +5726 @@
 
 !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) ::  edge   !<
-      CHARACTER(LEN=1), INTENT(IN) ::  var    !<
-
-      INTEGER(iwp) ::  i        !<
-      INTEGER(iwp) ::  ia       !<
-      INTEGER(iwp) ::  ib       !<
-      INTEGER(iwp) ::  ibgp     !<
-      INTEGER(iwp) ::  ijk      !<
-      INTEGER(iwp) ::  iw       !<
-      INTEGER(iwp) ::  j        !<
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 'l', 'r', 's' or 'n' 
+      CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
+
+      INTEGER(iwp) ::  i        !< Lower bound of the running index ia
+      INTEGER(iwp) ::  ia       !< i-index running over the parent-grid cell on the boundary 
+      INTEGER(iwp) ::  iaw      !< Reduced ia-index for workarr_lr
+      INTEGER(iwp) ::  iawbc    !< iaw-index pointing to the boundary-value nodes (either 0 or igsr-1)      
+      INTEGER(iwp) ::  ibc      !< Fixed i-index pointing to the boundary-value nodes (either i or iend)
+!AH      INTEGER(iwp) ::  ibeg     !< i-index pointing to the starting point of workarr_lr in the i-direction 
+      INTEGER(iwp) ::  iend     !< Upper bound of the running index ia
+      INTEGER(iwp) ::  ierr     !< MPI error code
+      INTEGER(iwp) ::  ijk      !< Running index for all child-grid cells within the anterpolation cell
+      INTEGER(iwp) ::  iw       !< i-index for wall_flags_0
+      INTEGER(iwp) ::  j        !< Running index in the y-direction
       INTEGER(iwp) ::  jco      !<
       INTEGER(iwp) ::  jcorr    !<
       INTEGER(iwp) ::  jinc     !<
-      INTEGER(iwp) ::  jw       !<
+      INTEGER(iwp) ::  jw       !< j-index for wall_flags_0
       INTEGER(iwp) ::  j1       !<
-      INTEGER(iwp) ::  k        !<
-      INTEGER(iwp) ::  k_wall   !< vertical index of topography top
+      INTEGER(iwp) ::  k        !< Running index in the z-direction
+      INTEGER(iwp) ::  k_wall   !< Vertical index of topography top
       INTEGER(iwp) ::  kco      !<
       INTEGER(iwp) ::  kcorr    !<
+      INTEGER(iwp) ::  kw       !< k-index for wall_flags_0
       INTEGER(iwp) ::  k1       !<
-      INTEGER(iwp) ::  l        !<
-      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        !<
-      REAL(wp) ::  fkpj        !<
-      REAL(wp) ::  fkpjp       !<
-      REAL(wp) ::  fk          !<
-      REAL(wp) ::  fkp         !<
-      REAL(wp) ::  rcorr       !<
-      REAL(wp) ::  rcorr_ijk   !<
- 
+      INTEGER(iwp) ::  l        !< Parent-grid running index in the x-direction
+      INTEGER(iwp) ::  lbeg     !< l-index pointing to the starting point of workarrc_lr in the l-direction 
+      INTEGER(iwp) ::  lp1      !< l+1
+      INTEGER(iwp) ::  loff     !< l-offset needed on the right boundary to correctly refer to boundary ghost points
+      INTEGER(iwp) ::  lw       !< Reduced l-index for workarrc_lr
+      INTEGER(iwp) ::  m        !< Parent-grid running index in the y-direction
+      INTEGER(iwp) ::  mnorthv  !< Upshift by one for the upper bound of index m in case of var == 'v'
+      INTEGER(iwp) ::  mp1      !< m+1
+      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  np1      !< n+1
+      INTEGER(iwp) ::  ntopw    !< Upshift by one for the upper bound of index n in case of var == 'w'
+      INTEGER(iwp) ::  var_flag !< Variable flag for BTEST( wall_flags_0 )
+
+      REAL(wp) ::  cellsum      !< Sum of child-grid node values over the anterpolation cell 
+      REAL(wP) ::  cellsumd     !< Sum of differences over the anterpolation cell 
+      REAL(wp) ::  fkj          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkjp         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpj         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpjp        !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fk           !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkp          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
+      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+
 ! 
 !--   Check which edge is to be handled
@@ -5534 +5786 @@
 !--      For u, nxl is a ghost node, but not for the other variables
          IF ( var == 'u' )  THEN
-            i  = nxl
-            ib = nxl - 1 
+            i     = nxl
+            iend  = nxl
+            ibc   = nxl
+            iawbc = 0
+            l     = icl + 2
+            lw    = 2
+            lbeg  = icl
+            loff  = 0
          ELSE
-            i  = nxl - 1
-            ib = nxl - 2
+            i     = nxl - igsr
+            iend  = nxl - 1
+            ibc   = nxl - 1
+            iawbc = igsr-1
+            l     = icl + 1
+            lw    = 1
+            lbeg  = icl
+            loff  = 0
          ENDIF
       ELSEIF ( edge == 'r' )  THEN
-         i  = nxr + 1
-         ib = nxr + 2
+         IF ( var == 'u' )  THEN
+            i     = nxr + 1            
+            iend  = nxr + 1
+            ibc   = nxr + 1 
+            iawbc = 0
+            l     = icr - 1
+            lw    = 1
+            lbeg  = icr - 2
+            loff  = 0
+         ELSE
+            i     = nxr + 1
+            iend  = nxr + igsr
+            ibc   = nxr + 1
+            iawbc = 0
+            l     = icr - 1
+            lw    = 1
+            lbeg  = icr - 2
+            loff  = 1
+         ENDIF
+      ENDIF
+
+      IF  ( var == 'w' )  THEN
+         ntopw = 1
+      ELSE
+         ntopw = 0
+      ENDIF
+
+      IF  ( var == 'v' )  THEN
+         mnorthv = 0
+      ELSE
+         mnorthv = 1
       ENDIF
 
@@ -5554 +5847 @@
          var_flag = 0
       ENDIF
-   
-      DO  j = nys, nyn+1
-!
-!--      Determine vertical index of topography top at grid point (j,i)
+
+      IF  ( var == 'u' )  THEN
+!AH!
+!AH!--      Substitute the necessary parent-grid data to the work array workarrc_lr. 
+!AH         workarrc_lr = 0.0_wp     
+!AH         IF  ( pdims(2) > 1 )  THEN
+!AH#if defined( __parallel )
+!AH            IF  ( nys == 0 )  THEN 
+!AH               workarrc_lr(0:cg%nz+1,jcsw:jcnw-1,0:2)                           &
+!AH                    = fc(0:cg%nz+1,jcsw:jcnw-1,lbeg:lbeg+2)
+!AH            ELSE IF  ( nyn == ny )  THEN
+!AH               workarrc_lr(0:cg%nz+1,jcsw+1:jcnw,0:2)                           &
+!AH                    = fc(0:cg%nz+1,jcsw+1:jcnw,lbeg:lbeg+2)
+!AH            ELSE
+!AH               workarrc_lr(0:cg%nz+1,jcsw+1:jcnw-1,0:2)                         &
+!AH                    = fc(0:cg%nz+1,jcsw+1:jcnw-1,lbeg:lbeg+2)
+!AH            ENDIF
+!AH!
+!AH!--         South-north exchange if more than one PE subdomain in the y-direction.
+!AH!--         Note that in case of 3-D nesting the south (psouth == MPI_PROC_NULL) 
+!AH!--         and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged 
+!AH!--         because the nest domain is not cyclic.
+!AH!--         From south to north
+!AH            CALL MPI_SENDRECV( workarrc_lr(0,jcsw+1,0), 1,                      &
+!AH                 workarrc_lr_exchange_type, psouth,  0,                         &
+!AH                 workarrc_lr(0,jcnw,0), 1,                                      &
+!AH                 workarrc_lr_exchange_type, pnorth,  0,                         &
+!AH                 comm2d, status, ierr ) 
+!AH!
+!AH!--         From north to south        
+!AH            CALL MPI_SENDRECV( workarrc_lr(0,jcnw-1,0), 1,                      &
+!AH                 workarrc_lr_exchange_type, pnorth,  1,                         &
+!AH                 workarrc_lr(0,jcsw,0), 1,                                      &
+!AH                 workarrc_lr_exchange_type, psouth,  1,                         &
+!AH                 comm2d, status, ierr )
+!AH#endif
+!AH         ELSE
+!AH            workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2)                                &
+!AH                 = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2)            
+!AH         ENDIF
+
+         DO  m = jcsw+1, jcnw-1
+            DO n = 0, kct 
+
+               DO  j = jfl(m), jfu(m)
+                  DO  k = kfl(n), kfu(n)
+                     f(k,j,ibc) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'v' )  THEN
+         
+         DO  m = jcsw+1, jcnw-1
+            DO n = 0, kct 
+
+               DO  j = jfl(m), jfl(m+1)-1
+                  DO  k = kfl(n), kfu(n) 
+                     f(k,j,ibc) = fc(n,m,l) 
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'w' )  THEN
+
+         DO  m = jcsw+1, jcnw-1
+            DO n = 1, kct + 1   ! It is important to go up to kct+1  
+
+               DO  j = jfl(m), jfu(m)
+                  f(nzb,j,ibc) = 0.0_wp   ! Because the n-loop starts from n=1 instead of 0
+                  DO  k = kfu(n-1)+1, kfu(n) 
+                     f(k,j,ibc) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE   ! scalars
+         
+         DO  m = jcsw+1, jcnw-1
+            DO n = 0, kct 
+               
+               DO  j = jfl(m), jfu(m)
+                  DO  k = kfl(n), kfu(n)
+                     f(k,j,ibc) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ENDIF  ! var
+
+   END SUBROUTINE pmci_interp_1sto_lr
+
+
+
+   SUBROUTINE pmci_interp_1sto_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
+                                   r2z, logc, logc_ratio, logc_kbounds,        &
+                                   nzt_topo_nestbc,                            &
+                                   kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc,   &
+                                   edge, var )
+
+!
+!--   Interpolation of ghost-node values used as the child-domain boundary
+!--   conditions. This subroutine handles the south and north boundaries. 
+!--   This subroutine is based on trilinear interpolation.
+
+      IMPLICIT NONE
+
+      INTEGER(iwp) ::  nzt_topo_nestbc   !<
+
+      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(1:2,0:ncorr-1,nzb:nzt_topo_nestbc,nxl:nxr),          &
+                                      INTENT(IN)    ::  logc_ratio    !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x           !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x           !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y           !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y           !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z           !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z           !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z           !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z           !<
+
+      
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)         ::  ic    !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)         ::  jc    !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), INTENT(IN)        ::  kc    !<
+      INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nxl:nxr),                &
+                                          INTENT(IN)           ::  logc  !<
+      INTEGER(iwp), DIMENSION(1:2,nxl:nxr), INTENT(IN)         ::  logc_kbounds  !< 
+
+      INTEGER(iwp) :: kct
+!AH
+!      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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH
+
+!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
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 'l', 'r', 's' or 'n' 
+      CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
+      
+      INTEGER(iwp) ::  i        !< Running index in the x-direction
+      INTEGER(iwp) ::  iinc     !<
+      INTEGER(iwp) ::  icorr    !<
+      INTEGER(iwp) ::  ico      !<
+      INTEGER(iwp) ::  ierr     !< MPI error code
+      INTEGER(iwp) ::  ijk      !< Running index for all child-grid cells within the anterpolation cell
+      INTEGER(iwp) ::  iw       !< i-index for wall_flags_0
+      INTEGER(iwp) ::  i1       !<
+      INTEGER(iwp) ::  j        !< Lower bound of the running index ja
+      INTEGER(iwp) ::  ja       !< Index in y-direction running over the parent-grid cell on the boundary 
+      INTEGER(iwp) ::  jaw      !< Reduced ja-index for workarr_sn
+      INTEGER(iwp) ::  jawbc    !< jaw-index pointing to the boundary-value nodes (either 0 or jgsr-1)
+!AH      INTEGER(iwp) ::  jbeg     !< j-index pointing to the starting point of workarr_sn in the j-direction
+      INTEGER(iwp) ::  jbc      !< Fixed j-index pointing to the boundary-value nodes (either j or jend)
+      INTEGER(iwp) ::  jend     !< Upper bound of the running index ja
+      INTEGER(iwp) ::  jw       !< j-index for wall_flags_0
+      INTEGER(iwp) ::  k        !< Running index in the z-direction
+      INTEGER(iwp) ::  k_wall   !< Vertical index of topography top
+      INTEGER(iwp) ::  kcorr    !<
+      INTEGER(iwp) ::  kco      !<
+      INTEGER(iwp) ::  kw       !< k-index for wall_flags_0
+      INTEGER(iwp) ::  k1       !<
+      INTEGER(iwp) ::  l        !< Parent-grid running index in the x-direction
+      INTEGER(iwp) ::  lp1      !< l+1
+      INTEGER(iwp) ::  lrightu  !< Upshift by one for the upper bound of index l in case of var == 'u'
+      INTEGER(iwp) ::  m        !< Parent-grid running index in the y-direction
+      INTEGER(iwp) ::  mbeg     !< m-index pointing to the starting point of workarrc_sn in the m-direction 
+      INTEGER(iwp) ::  moff     !< m-offset needed on the north boundary to correctly refer to boundary ghost points
+      INTEGER(iwp) ::  mp1      !< m+1
+      INTEGER(iwp) ::  mw       !< Reduced m-index for workarrc_sn
+      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  np1      !< n+1
+      INTEGER(iwp) ::  ntopw    !< Upshift by one for the upper bound of index n in case of var == 'w'
+      INTEGER(iwp) ::  var_flag !< Variable flag for BTEST( wall_flags_0 )
+
+      REAL(wp) ::  cellsum      !< Sum of child-grid node values over the anterpolation cell 
+      REAL(wp) ::  cellsumd     !< Sum of differences over the anterpolation cell 
+      REAL(wp) ::  fk           !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkj          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkjp         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpj         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpjp        !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkp          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
+      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+
+!
+!--   Check which edge is to be handled: south or north
+      IF ( edge == 's' )  THEN
+!
+!--      For v, nys is a ghost node, but not for the other variables
+         IF ( var == 'v' )  THEN
+            j     = nys
+            jend  = nys 
+            jawbc = 0
+            jbc   = nys
+            m     = jcs + 2
+            mw    = 2
+            mbeg  = jcs
+            moff  = 0
+         ELSE
+            j     = nys - jgsr
+            jend  = nys - 1
+            jawbc = jgsr - 1
+            jbc   = nys - 1
+            m     = jcs + 1
+            mw    = 1
+            mbeg  = jcs
+            moff  = 0
+         ENDIF
+      ELSEIF ( edge == 'n' )  THEN
+         IF ( var == 'v' )  THEN
+            j     = nyn + 1
+            jend  = nyn + 1
+            jawbc = 0
+            jbc   = nyn + 1
+            m     = jcn - 1
+            mw    = 1
+            mbeg  = jcn - 2
+            moff  = 0
+         ELSE
+            j     = nyn + 1
+            jend  = nyn + jgsr
+            jawbc = 0
+            jbc   = nyn + 1
+            m     = jcn - 1
+            mw    = 1
+            mbeg  = jcn - 2
+            moff  = 1
+         ENDIF
+      ENDIF
+
+      IF  ( var == 'w' )  THEN
+         ntopw = 1
+      ELSE
+         ntopw = 0
+      ENDIF
+
+      IF  ( var == 'u' )  THEN
+         lrightu = 0
+      ELSE
+         lrightu = 1
+      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
+
+      IF  ( var == 'v' )  THEN
+!AH!
+!AH!--      Substitute the necessary parent-grid data to the work array workarrc_sn. 
+!AH         workarrc_sn = 0.0_wp     
+!AH         IF  ( pdims(1) > 1 )  THEN
+!AH#if defined( __parallel )
+!AH            IF  ( nxl == 0 )  THEN   ! if ( bc_dirichlet_l )
+!AH               workarrc_sn(0:cg%nz+1,0:2,iclw:icrw-1)                           &
+!AH                    = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1)
+!AH            ELSE IF  ( nxr == nx )  THEN    ! if ( bc_dirichlet_r )
+!AH               workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw)                           &
+!AH                    = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw)
+!AH            ELSE
+!AH               workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                         &
+!AH                    = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
+!AH            ENDIF
+!AH!
+!AH!--         Left-right exchange if more than one PE subdomain in the x-direction.
+!AH!--         Note that in case of 3-D nesting the left (pleft == MPI_PROC_NULL) and 
+!AH!--         right (pright == MPI_PROC_NULL) boundaries are not exchanged because 
+!AH!--         the nest domain is not cyclic.
+!AH!--         From left to right
+!AH            CALL MPI_SENDRECV( workarrc_sn(0,0,iclw+1), 1,                      &
+!AH                 workarrc_sn_exchange_type, pleft,   0,                         &
+!AH                 workarrc_sn(0,0,icrw), 1,                                      &
+!AH                 workarrc_sn_exchange_type, pright,  0,                         &
+!AH                 comm2d, status, ierr ) 
+!AH!
+!AH!--         From right to left        
+!AH            CALL MPI_SENDRECV( workarrc_sn(0,0,icrw-1), 1,                      &
+!AH                 workarrc_sn_exchange_type, pright,  1,                         &
+!AH                 workarrc_sn(0,0,iclw), 1,                                      &
+!AH                 workarrc_sn_exchange_type, pleft,   1,                         &
+!AH                 comm2d, status, ierr )
+!AH#endif      
+!AH         ELSE
+!AH            workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                            &
+!AH                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
+!AH         ENDIF
+
+         DO  l = iclw+1, icrw-1
+            DO n = 0, kct 
+
+               DO  i = ifl(l), ifu(l)
+                  DO  k = kfl(n), kfu(n)
+                     f(k,jbc,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'u' )  THEN
+         
+         DO  l = iclw+1, icrw-1
+            DO n = 0, kct 
+               
+               DO  i = ifl(l), ifl(l+1)-1
+                  DO  k = kfl(n), kfu(n) 
+                     f(k,jbc,i) = fc(n,m,l) 
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'w' )  THEN
+
+         DO  l = iclw+1, icrw-1
+            DO n = 1, kct + 1   ! It is important to go up to kct+1  
+               
+               DO  i = ifl(l), ifu(l)
+                  f(nzb,jbc,i) = 0.0_wp   ! Because the n-loop starts from n=1 instead of 0
+                  DO  k = kfu(n-1)+1, kfu(n) 
+                     f(k,jbc,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE   ! scalars
+         
+         DO  l = iclw+1, icrw-1
+            DO n = 0, kct 
+               
+               DO  i = ifl(l), ifu(l)
+                  DO  k = kfl(n), kfu(n)
+                     f(k,jbc,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ENDIF  ! var
+
+   END SUBROUTINE pmci_interp_1sto_sn
+
+
+
+   SUBROUTINE pmci_interp_1sto_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,       &
+                                  r1z, r2z, kct, ifl, ifu, jfl, jfu, kfl, kfu, &
+                                  ijkfc, var )
+
+!
+!--   Interpolation of ghost-node values used as the child-domain boundary
+!--   conditions. This subroutine handles the top boundary. 
+!--   This subroutine is based on trilinear interpolation.
+
+      IMPLICIT NONE
+
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
+                                      INTENT(INOUT) ::  f     !< Child-grid array
+      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr),                          &
+                                      INTENT(IN)    ::  fc    !< Parent-grid array
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x   !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x   !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y   !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y   !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z   !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z   !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z   !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z   !<
+
+      
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  ic    !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  jc    !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  kc    !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  kc    !<
+
+      INTEGER(iwp) :: kct
+!AH
+!      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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), 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
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) :: var   !<
+
+      INTEGER(iwp) ::  i    !<
+      INTEGER(iwp) ::  ib   !<
+      INTEGER(iwp) ::  iclc !< Lower i-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  icrc !< Upper i-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  ie   !<
+      INTEGER(iwp) ::  ierr !< MPI error code
+      INTEGER(iwp) ::  ijk  !<
+      INTEGER(iwp) ::  iw   !<
+      INTEGER(iwp) ::  j    !<
+      INTEGER(iwp) ::  jb   !<
+      INTEGER(iwp) ::  jcsc !< Lower j-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  jcnc !< Upper j-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  je   !< 
+      INTEGER(iwp) ::  jw   !<     
+      INTEGER(iwp) ::  k    !< Vertical child-grid index fixed to the boundary-value level 
+      INTEGER(iwp) ::  ka   !< Running vertical child-grid index
+      INTEGER(iwp) ::  kw   !< 
+      INTEGER(iwp) ::  l    !< Parent-grid index in x-direction
+      INTEGER(iwp) ::  lp1  !< l+1
+      INTEGER(iwp) ::  m    !< Parent-grid index in y-direction
+      INTEGER(iwp) ::  mp1  !< m+1
+      INTEGER(iwp) ::  n    !< Parent-grid work array index in z-direction
+      INTEGER(iwp) ::  np1  !< n+1
+      INTEGER(iwp) ::  noff !< n-offset needed on the top boundary to correctly refer to boundary ghost points  
+      INTEGER(iwp) ::  nw   !< n-index for workarrc_t
+      INTEGER(iwp) ::  var_flag  !<
+      
+      REAL(wp) ::  cellsum     !<
+      REAL(wp) ::  cellsumd    !<
+      REAL(wp) ::  fac
+      REAL(wp) ::  fk          !<
+      REAL(wp) ::  fkj         !<
+      REAL(wp) ::  fkjp        !<
+      REAL(wp) ::  fkpj        !<
+      REAL(wp) ::  fkpjp       !<
+      REAL(wp) ::  fkp         !<
+      REAL(wp) ::  rcorr       !<
+      REAL(wp) ::  rcorr_ijk   !<
+
+
+      IF ( var == 'w' )  THEN
+         k    = nzt
+         noff = 0
+      ELSE
+         k    = nzt + 1
+         noff = 1
+      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
+      n  = kc(k) + noff
+      nw = 1 
+
+      IF  ( var == 'w' )  THEN
+!AH!
+!AH!--      Substitute the necessary parent-grid data to the work array. 
+!AH!--      Note that the dimension of workarrc_t is (0:2,jcsw:jcnw,iclw:icrw),
+!AH!--      And the jc?w and ic?w-index bounds depend on the location of the PE-
+!AH!--      subdomain relative to the side boundaries. 
+!AH         iclc = iclw + 1
+!AH         icrc = icrw - 1      
+!AH         jcsc = jcsw + 1
+!AH         jcnc = jcnw - 1
+!AH         IF  ( bc_dirichlet_l )  THEN
+!AH            iclc = iclw
+!AH         ENDIF
+!AH         IF  ( bc_dirichlet_r )  THEN
+!AH            icrc = icrw
+!AH         ENDIF
+!AH         IF  ( bc_dirichlet_s )  THEN
+!AH            jcsc = jcsw
+!AH         ENDIF
+!AH         IF  ( bc_dirichlet_n )  THEN
+!AH            jcnc = jcnw
+!AH         ENDIF
+!AH         workarrc_t = 0.0_wp
+!AH         workarrc_t(0:2,jcsc:jcnc,iclc:icrc)                                    &
+!AH              = fc(kc(k):kc(k)+2,jcsc:jcnc,iclc:icrc)
+!AH!
+!AH!--      Left-right exchange if more than one PE subdomain in the x-direction.
+!AH!--      Note that in case of 3-D nesting the left and right boundaries are 
+!AH!--      not exchanged because the nest domain is not cyclic.
+!AH#if defined( __parallel )
+!AH         IF  ( pdims(1) > 1 )  THEN
+!AH!
+!AH!--         From left to right
+!AH            CALL MPI_SENDRECV( workarrc_t(0,jcsw,iclw+1), 1,                    &
+!AH                 workarrc_t_exchange_type_y, pleft,  0,                         &
+!AH                 workarrc_t(0,jcsw,icrw), 1,                                    &
+!AH                 workarrc_t_exchange_type_y, pright, 0,                         &
+!AH                 comm2d, status, ierr ) 
+!AH!
+!AH!--         From right to left        
+!AH            CALL MPI_SENDRECV( workarrc_t(0,jcsw,icrw-1), 1,                    &
+!AH                 workarrc_t_exchange_type_y, pright, 1,                         &
+!AH                 workarrc_t(0,jcsw,iclw), 1,                                    &
+!AH                 workarrc_t_exchange_type_y, pleft,  1,                         &
+!AH                 comm2d, status, ierr )
+!AH         ENDIF
+!AH!
+!AH!--      South-north exchange if more than one PE subdomain in the y-direction.
+!AH!--      Note that in case of 3-D nesting the south and north boundaries are 
+!AH!--      not exchanged because the nest domain is not cyclic.
+!AH         IF  ( pdims(2) > 1 )  THEN
+!AH!
+!AH!--         From south to north         
+!AH            CALL MPI_SENDRECV( workarrc_t(0,jcsw+1,iclw), 1,                    &
+!AH                 workarrc_t_exchange_type_x, psouth, 2,                         &
+!AH                 workarrc_t(0,jcnw,iclw), 1,                                    &
+!AH                 workarrc_t_exchange_type_x, pnorth, 2,                         &
+!AH                 comm2d, status, ierr ) 
+!AH!
+!AH!--         From north to south        
+!AH            CALL MPI_SENDRECV( workarrc_t(0,jcnw-1,iclw), 1,                    &
+!AH                 workarrc_t_exchange_type_x, pnorth, 3,                         &
+!AH                 workarrc_t(0,jcsw,iclw), 1,                                    &
+!AH                 workarrc_t_exchange_type_x, psouth, 3,                         &
+!AH                 comm2d, status, ierr )
+!AH         ENDIF
+!AH#endif      
+!AH
+
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     f(k,j,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'u' )  THEN
+         
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+               
+               DO  i = ifl(l), ifl(l+1)-1
+                  DO  j = jfl(m), jfu(m) 
+                     f(k,j,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE IF  ( var == 'v' )  THEN
+
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfl(m+1)-1
+                     f(k,j,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ELSE   ! scalars
+         
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     f(k,j,i) = fc(n,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+      ENDIF  ! var
+
+   END SUBROUTINE pmci_interp_1sto_t
+
+
+
+   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
+                                   r2z, logc, logc_ratio, logc_kbounds,        &
+                                   nzt_topo_nestbc,                            &
+                                   kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc,   &
+                                   edge, var )
+!
+!--   Interpolation of ghost-node values used as the child-domain boundary
+!--   conditions. This subroutine handles the left and right boundaries. It is
+!--   based on trilinear interpolation.
+
+      IMPLICIT NONE
+
+      INTEGER(iwp) ::  nzt_topo_nestbc   !<
+
+      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(1:2,0:ncorr-1,nzb:nzt_topo_nestbc,nys:nyn),          &
+                                      INTENT(IN)    ::  logc_ratio   !<
+!AH      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r1x     !<
+!AH      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r2x     !<
+!AH      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r1y     !<
+!AH      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r2y     !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x     !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x     !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y     !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y     !<
+
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z     !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z     !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z     !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z     !<
+      
+
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)         ::  ic     !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)         ::  jc     !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), 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 !<
+
+      INTEGER(iwp) :: kct
+
+!AH
+!      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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH
+
+!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
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 'l', 'r', 's' or 'n' 
+      CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
+
+      INTEGER(iwp) ::  i        !< Lower bound of the running index ia
+      INTEGER(iwp) ::  ia       !< i-index running over the parent-grid cell on the boundary 
+      INTEGER(iwp) ::  iaw      !< Reduced ia-index for workarr_lr
+      INTEGER(iwp) ::  iawbc    !< iaw-index pointing to the boundary-value nodes (either 0 or igsr-1)      
+      INTEGER(iwp) ::  ibc      !< Fixed i-index pointing to the boundary-value nodes (either i or iend)
+!AH      INTEGER(iwp) ::  ibeg     !< i-index pointing to the starting point of workarr_lr in the i-direction 
+      INTEGER(iwp) ::  iend     !< Upper bound of the running index ia
+      INTEGER(iwp) ::  ierr     !< MPI error code
+      INTEGER(iwp) ::  ijk      !< Running index for all child-grid cells within the anterpolation cell
+      INTEGER(iwp) ::  iw       !< i-index for wall_flags_0
+      INTEGER(iwp) ::  j        !< Running index in the y-direction
+      INTEGER(iwp) ::  jco      !<
+      INTEGER(iwp) ::  jcorr    !<
+      INTEGER(iwp) ::  jinc     !<
+      INTEGER(iwp) ::  jw       !< j-index for wall_flags_0
+      INTEGER(iwp) ::  j1       !<
+      INTEGER(iwp) ::  k        !< Running index in the z-direction
+      INTEGER(iwp) ::  k_wall   !< Vertical index of topography top
+      INTEGER(iwp) ::  kco      !<
+      INTEGER(iwp) ::  kcorr    !<
+      INTEGER(iwp) ::  kw       !< k-index for wall_flags_0
+      INTEGER(iwp) ::  k1       !<
+      INTEGER(iwp) ::  l        !< Parent-grid running index in the x-direction
+      INTEGER(iwp) ::  lbeg     !< l-index pointing to the starting point of workarrc_lr in the l-direction 
+      INTEGER(iwp) ::  lp1      !< l+1
+      INTEGER(iwp) ::  loff     !< l-offset needed on the right boundary to correctly refer to boundary ghost points
+      INTEGER(iwp) ::  lw       !< Reduced l-index for workarrc_lr
+      INTEGER(iwp) ::  m        !< Parent-grid running index in the y-direction
+      INTEGER(iwp) ::  mnorthv  !< Upshift by one for the upper bound of index m in case of var == 'v'
+      INTEGER(iwp) ::  mp1      !< m+1
+      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  np1      !< n+1
+      INTEGER(iwp) ::  ntopw    !< Upshift by one for the upper bound of index n in case of var == 'w'
+      INTEGER(iwp) ::  var_flag !< Variable flag for BTEST( wall_flags_0 )
+
+      REAL(wp) ::  cellsum      !< Sum of child-grid node values over the anterpolation cell 
+      REAL(wP) ::  cellsumd     !< Sum of differences over the anterpolation cell 
+      REAL(wp) ::  fkj          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkjp         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpj         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpjp        !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fk           !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkp          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
+      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+ 
+! 
+!--   Check which edge is to be handled
+      IF ( edge == 'l' )  THEN
+!
+!--      For u, nxl is a ghost node, but not for the other variables
+         IF ( var == 'u' )  THEN
+            i     = nxl
+            iend  = nxl
+            ibc   = nxl
+            iawbc = 0
+            lbeg  = icl
+            loff  = 0
+         ELSE
+            i     = nxl - igsr
+            iend  = nxl - 1
+            ibc   = nxl - 1
+            iawbc = igsr-1
+            lbeg  = icl
+            loff  = 0
+         ENDIF
+      ELSEIF ( edge == 'r' )  THEN
+         IF ( var == 'u' )  THEN
+            i     = nxr + 1            
+            iend  = nxr + 1
+            ibc   = nxr + 1 
+            iawbc = 0
+            lbeg  = icr - 2
+            loff  = 0
+         ELSE
+            i     = nxr + 1
+            iend  = nxr + igsr
+            ibc   = nxr + 1
+            iawbc = 0
+            lbeg  = icr - 2
+            loff  = 1
+         ENDIF
+
+      ENDIF
+
+      IF  ( var == 'w' )  THEN
+         ntopw = 1
+      ELSE
+         ntopw = 0
+      ENDIF
+
+      IF  ( var == 'v' )  THEN
+         mnorthv = 0
+      ELSE
+         mnorthv = 1
+      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
+
+!AH
+!
+!--   Substitute the necessary parent-grid data to the work array workarrc_lr. 
+      workarrc_lr = 0.0_wp     
+      IF  ( pdims(2) > 1 )  THEN
+#if defined( __parallel )
+         IF  ( nys == 0 )  THEN 
+            workarrc_lr(0:cg%nz+1,jcsw:jcnw-1,0:2)                              &
+                 = fc(0:cg%nz+1,jcsw:jcnw-1,lbeg:lbeg+2)
+         ELSE IF  ( nyn == ny )  THEN
+            workarrc_lr(0:cg%nz+1,jcsw+1:jcnw,0:2)                              &
+                 = fc(0:cg%nz+1,jcsw+1:jcnw,lbeg:lbeg+2)
+         ELSE
+            workarrc_lr(0:cg%nz+1,jcsw+1:jcnw-1,0:2)                            &
+                 = fc(0:cg%nz+1,jcsw+1:jcnw-1,lbeg:lbeg+2)
+         ENDIF
+!
+!--      South-north exchange if more than one PE subdomain in the y-direction.
+!--      Note that in case of 3-D nesting the south (psouth == MPI_PROC_NULL) 
+!--      and north (pnorth == MPI_PROC_NULL) boundaries are not exchanged 
+!--      because the nest domain is not cyclic.
+!--      From south to north
+         CALL MPI_SENDRECV( workarrc_lr(0,jcsw+1,0), 1,                         &
+              workarrc_lr_exchange_type, psouth,  0,                            &
+              workarrc_lr(0,jcnw,0), 1,                                         &
+              workarrc_lr_exchange_type, pnorth,  0,                            &
+              comm2d, status, ierr ) 
+!
+!--      From north to south        
+         CALL MPI_SENDRECV( workarrc_lr(0,jcnw-1,0), 1,                         &
+              workarrc_lr_exchange_type, pnorth,  1,                            &
+              workarrc_lr(0,jcsw,0), 1,                                         &
+              workarrc_lr_exchange_type, psouth,  1,                            &
+              comm2d, status, ierr )
+#endif
+      ELSE
+         workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2)                                   &
+              = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2)            
+      ENDIF
+!
+!AH
+      workarr_lr = 0.0_wp
+      
+      DO  ia = i, iend
+         iaw = ia - i
+         DO  j = nys-1, nyn+1
+!AH         DO  j = nys, nyn
+!
+!--         Determine vertical index of topography top at grid point (j,i)
 !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)
-            n = kc(k)
-            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
-            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
-            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
-            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
-            fk       = r1y(j) * fkj  + r2y(j) * fkjp
-            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
-            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+            DO  k = nzb, nzt+1 !k_wall, nzt+1
+!AH               l   = ic(ia) - ic(i)
+               l   = ic(ia) - lbeg
+               lp1 = MIN( l + 1, 2 )  ! If l+1 > 2 (l=ic(nxr+1)-lbeg), r1x = 1 and r2x = 0
+               m   = jc(j)
+               mp1 = MIN( m + 1, jcnw )  ! If m+1 > jcn (m=jc(nyn+1)), r1y = 1 and r2y = 0
+               n   = kc(k)
+               np1 = n + 1            
+
+!AH
+!               fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
+!               fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
+!               fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
+!               fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
+!AH
+               
+               fkj      = r1x(ia) * workarrc_lr(n,m,l)     + r2x(ia) * workarrc_lr(n,m,lp1)
+               fkjp     = r1x(ia) * workarrc_lr(n,mp1,l)   + r2x(ia) * workarrc_lr(n,mp1,lp1)
+               fkpj     = r1x(ia) * workarrc_lr(np1,m,l)   + r2x(ia) * workarrc_lr(np1,m,lp1)
+               fkpjp    = r1x(ia) * workarrc_lr(np1,mp1,l) + r2x(ia) * workarrc_lr(np1,mp1,lp1)
+               
+               fk       = r1y(j) * fkj  + r2y(j) * fkjp
+               fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
+
+!AH
+!               f(k,j,i) = r1z(k) * fk + r2z(k) * fkp               
+               workarr_lr(k,j,iaw) = r1z(k) * fk + r2z(k) * fkp
+               
+!               if  ( ( edge == 'l' ) .and. ( ia  == ibc ) ) then
+!                  write(9,"('pmci_interp_tril_lr: ',a2,2x,11(i4,2x),7(e12.5,2x))") var, k, j, ia, ibc, iaw, n, m, l, np1, mp1, lp1, &
+!                       workarr_lr(k,j,iaw), r1x(ia), r2x(ia), workarrc_lr(n,m,l), workarrc_lr(n,mp1,l), workarrc_lr(np1,m,l), workarrc_lr(np1,mp1,l)
+!                  flush(9)
+!               endif
+
+!AH
+               
+            ENDDO
          ENDDO
       ENDDO
@@ -5651 +6813 @@
       ENDIF  ! ( topography /= 'flat' )
 !
-!--   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
+!--   Apply the reversibility correction. 
+
+!      if  ( var == 'u' )  then
+
+      l  = ic(ibc) + loff
+      lw = 1
+!      write(9,"('pmci_interp_tril_lr: edge, var, l, i, iend, ifl(l), ifu(l) = ',2(a2,2x),5(i4,2x))") &
+!           edge, var, l, i, iend, ifl(l), ifu(l)
+!      flush(9)
+!AH      DO  m = jcsw+1, jcnw-1       
+      DO  m = jcsw + 1, jcnw - mnorthv   ! mnorthv = 0 for v and 1 for all others
+         DO  n = 0, kct + ntopw            ! ntopw = 1 for w and 0 for all others
+            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)
+               iaw = ia - i
+               iw  = MAX( MIN( ia, nx+1 ), -1 )
+               DO  j = jfl(m), jfu(m)
+                  jw = MAX( MIN( j, ny+1 ), -1 )
+                  DO  k = kfl(n), kfu(n)
+                     kw = MIN( k, nzt+1 )
+!AH
+!                     cellsum = cellsum + MERGE( f(k,j,ia), 0.0_wp,              &
+!                          BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+                     cellsum = cellsum + MERGE( workarr_lr(k,j,iaw), 0.0_wp,     &
+                          BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+!AH
+                     
+!AH                     celltmpd(ijk) = ABS( fc(n,m,l) - f(k,j,ia) )
+!AH                     celltmpd(ijk) = ABS( workarrc_lr(n,m,lw) - f(k,j,ia) )
+                     celltmpd(ijk) = ABS( workarrc_lr(n,m,lw) - workarr_lr(k,j,iaw) )
+                     cellsumd      = cellsumd  + MERGE( celltmpd(ijk),          &
+                          0.0_wp, BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+
+!                     write(9,"('lr1: ',a1,2x,8(i4,2x),5(e12.5,2x))") var, n, m, l, k, j, ia, iaw, ijk,  &
+!                          workarrc_lr(n,m,lw), workarr_lr(k,j,iaw), cellsum, celltmpd(ijk), cellsumd
+!                     flush(9)
+
+                     ijk = ijk + 1
                   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
+            ENDDO
+
+            IF ( ijkfc(n,m,l) /= 0 )  THEN
+               cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+!               rcorr     = fc(n,m,l) - cellsum
+               rcorr     = workarrc_lr(n,m,lw) - 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
+!--         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)
+!AH               iaw = ia - ifl(l)
+               iaw = ia - i
+               DO  j = jfl(m), jfu(m)
+                  DO  k = kfl(n), kfu(n)
+                     rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+!AH                     f(k,j,ia) = f(k,j,ia) + rcorr_ijk
+                     workarr_lr(k,j,iaw) = workarr_lr(k,j,iaw) + rcorr_ijk
+
+!                     write(9,"('lr2: ',a1,2x,9(i4,2x),4(e12.5,2x))") var, n, m, l, k, j, ia, iaw, ijk, ijkfc(n,m,l), &
+!                          rcorr, rcorr_ijk, workarr_lr(k,j,iaw), workarrc_lr(n,m,lw)
+!                     flush(9)
+
+                     ijk = ijk + 1
                   ENDDO
                ENDDO
+            ENDDO
+!
+!--         Velocity components tangential to the boundary need special 
+!--         treatment because their anterpolation cells are flat in their 
+!--         direction and hence do not cover all the child-grid boundary nodes. 
+!--         These components are next linearly interpolated to those child-grid 
+!--         boundary nodes which are not covered by the anterpolation cells.
+            IF  ( ( var == 'w' ) .AND. ( n > 0 ) )  THEN
+               DO  k = kfu(n-1)+1, kfl(n)-1
+                  IF  ( k <= nzt )  THEN
+                     DO  j = jfl(m), jfu(m)
+                        IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+!AH                              
+!                           f(k,j,ia) = r1z(k) * f(kfu(n-1),j,ia)                &
+!                                   + r2z(k) * f(kfl(n),j,ia)
+                           workarr_lr(k,j,iawbc) =                              &
+                                  r1z(k) * workarr_lr(kfu(n-1),j,iawbc)         &
+                                + r2z(k) * workarr_lr(kfl(n),j,iawbc)
+
+!                           write(9,"('lr3: ',a1,2x,7(i4,2x),3(e12.5,2x))") var, n, m, l, k, j, iawbc, ibc, &
+!                                workarr_lr(k-1,j,iawbc), workarr_lr(k,j,iawbc), workarr_lr(k+1,j,iawbc)
+!                           flush(9)
+
+!AH                              
+                        ENDIF
+                     ENDDO
+                  ENDIF
+               ENDDO
+            ENDIF
+
+            IF  ( ( var == 'v' ) .AND. ( m > jcsw ) )  THEN
+               DO  j = jfu(m-1)+1, jfl(m)-1
+                  IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+                     DO  k = kfl(n), kfu(n)
+                        IF  ( k <= nzt+1 )  THEN                              
+!AH                              
+!                           f(k,j,ia) = r1y(j) * f(k,jfu(m-1),ia)                &
+!                                   + r2y(j) * f(k,jfl(m),ia)
+                           workarr_lr(k,j,iawbc) =                              &
+                                  r1y(j) * workarr_lr(k,jfu(m-1),iawbc)         &
+                                + r2y(j) * workarr_lr(k,jfl(m),iawbc)
+
+!                           write(9,"('lr4: ',a1,2x,7(i4,2x),3(e12.5,2x))") var, n, m, l, k, j, iawbc, ibc, &
+!                                workarr_lr(k,j-1,iawbc), workarr_lr(k,j,iawbc), workarr_lr(k,j+1,iawbc)
+!                           flush(9)
+
+!AH                              
+                        ENDIF
+                     ENDDO
+                  ENDIF
+               ENDDO
+            ENDIF
                
-            ENDDO  ! n
-         ENDDO  ! m
-         
-      ENDIF  ! var not v or w 
-!
-!--   Store the boundary values also into the other redundant ghost node layers.
-!--   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
-         DO  ibgp = -nbgp, ib
-            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
-         ENDDO
-      ELSEIF ( edge == 'r' )  THEN
-         DO  ibgp = ib, nx+nbgp
-            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
-         ENDDO
-      ENDIF
+         ENDDO  ! n
+      ENDDO  ! m
+
+!      endif  ! var
+
+!
+!--   Finally substitute the boundary values.
+      f(nzb:nzt+1,nys:nyn,ibc) = workarr_lr(nzb:nzt+1,nys:nyn,iawbc)
+
+!      do  k = 0, 2
+!         do  j = nys, nyn 
+!            if  ( edge == 'l' )  then
+!               write(9,"('lr5: ',2(a2,2x),4(i4,2x),4(e12.5,2x))") edge, var, k, j, ibc, iawbc,  &
+!                    workarr_lr(k,j,iawbc), f(k,j,ibc), f(k,j,ibc+1), f(k,j,ibc+2)
+!            else if  ( edge == 'r' )  then
+!               write(9,"('lr5: ',2(a2,2x),4(i4,2x),4(e12.5,2x))") edge, var, k, j, ibc, iawbc,  &
+!                    f(k,j,ibc-2), f(k,j,ibc-1), f(k,j,ibc), workarr_lr(k,j,iawbc)
+!            endif
+!         enddo
+!      enddo
+!      flush(9)
 
    END SUBROUTINE pmci_interp_tril_lr
@@ -5754 +7000 @@
       REAL(wp), DIMENSION(1:2,0:ncorr-1,nzb:nzt_topo_nestbc,nxl:nxr),          &
                                       INTENT(IN)    ::  logc_ratio    !<
-      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           !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x           !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x           !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y           !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y           !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z           !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z           !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z           !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z           !<
+
       
-      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !<
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !<
-      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)         ::  ic    !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)         ::  jc    !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), INTENT(IN)        ::  kc    !<
       INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nxl:nxr),                &
                                           INTENT(IN)           ::  logc  !<
@@ -5769 +7019 @@
 
       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(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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH
+
 !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
@@ -5778 +7035 @@
       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) ::  edge   !<
-      CHARACTER(LEN=1), INTENT(IN) ::  var    !<
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 'l', 'r', 's' or 'n' 
+      CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
       
-      INTEGER(iwp) ::  i       !<
-      INTEGER(iwp) ::  iinc    !<
-      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) ::  kcorr   !<
-      INTEGER(iwp) ::  kco     !<
-      INTEGER(iwp) ::  k1      !<
-      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         !<
-      REAL(wp) ::  fkjp        !<
-      REAL(wp) ::  fkpj        !<
-      REAL(wp) ::  fkpjp       !<
-      REAL(wp) ::  fkp         !<
-      REAL(wp) ::  rcorr       !< 
-      REAL(wp) ::  rcorr_ijk   !< 
-      
+      INTEGER(iwp) ::  i        !< Running index in the x-direction
+      INTEGER(iwp) ::  iinc     !<
+      INTEGER(iwp) ::  icorr    !<
+      INTEGER(iwp) ::  ico      !<
+      INTEGER(iwp) ::  ierr     !< MPI error code
+      INTEGER(iwp) ::  ijk      !< Running index for all child-grid cells within the anterpolation cell
+      INTEGER(iwp) ::  iw       !< i-index for wall_flags_0
+      INTEGER(iwp) ::  i1       !<
+      INTEGER(iwp) ::  j        !< Lower bound of the running index ja
+      INTEGER(iwp) ::  ja       !< Index in y-direction running over the parent-grid cell on the boundary 
+      INTEGER(iwp) ::  jaw      !< Reduced ja-index for workarr_sn
+      INTEGER(iwp) ::  jawbc    !< jaw-index pointing to the boundary-value nodes (either 0 or jgsr-1)
+!AH      INTEGER(iwp) ::  jbeg     !< j-index pointing to the starting point of workarr_sn in the j-direction
+      INTEGER(iwp) ::  jbc      !< Fixed j-index pointing to the boundary-value nodes (either j or jend)
+      INTEGER(iwp) ::  jend     !< Upper bound of the running index ja
+      INTEGER(iwp) ::  jw       !< j-index for wall_flags_0
+      INTEGER(iwp) ::  k        !< Running index in the z-direction
+      INTEGER(iwp) ::  k_wall   !< Vertical index of topography top
+      INTEGER(iwp) ::  kcorr    !<
+      INTEGER(iwp) ::  kco      !<
+      INTEGER(iwp) ::  kw       !< k-index for wall_flags_0
+      INTEGER(iwp) ::  k1       !<
+      INTEGER(iwp) ::  l        !< Parent-grid running index in the x-direction
+      INTEGER(iwp) ::  lp1      !< l+1
+      INTEGER(iwp) ::  lrightu  !< Upshift by one for the upper bound of index l in case of var == 'u'
+      INTEGER(iwp) ::  m        !< Parent-grid running index in the y-direction
+      INTEGER(iwp) ::  mbeg     !< m-index pointing to the starting point of workarrc_sn in the m-direction 
+      INTEGER(iwp) ::  moff     !< m-offset needed on the north boundary to correctly refer to boundary ghost points
+      INTEGER(iwp) ::  mp1      !< m+1
+      INTEGER(iwp) ::  mw       !< Reduced m-index for workarrc_sn
+      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  np1      !< n+1
+      INTEGER(iwp) ::  ntopw    !< Upshift by one for the upper bound of index n in case of var == 'w'
+      INTEGER(iwp) ::  var_flag !< Variable flag for BTEST( wall_flags_0 )
+
+      REAL(wp) ::  cellsum      !< Sum of child-grid node values over the anterpolation cell 
+      REAL(wp) ::  cellsumd     !< Sum of differences over the anterpolation cell 
+      REAL(wp) ::  fk           !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkj          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkjp         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpj         !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkpjp        !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  fkp          !< Intermediate result in trilinear interpolation
+      REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
+      REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+
 !
 !--   Check which edge is to be handled: south or north
@@ -5820 +7095 @@
 !--      For v, nys is a ghost node, but not for the other variables
          IF ( var == 'v' )  THEN
-            j  = nys
-            jb = nys - 1 
+            j     = nys
+            jend  = nys 
+            jawbc = 0
+            jbc   = nys       
+            mbeg  = jcs
+            moff  = 0
          ELSE
-            j  = nys - 1
-            jb = nys - 2
+            j     = nys - jgsr
+            jend  = nys - 1
+            jawbc = jgsr - 1
+            jbc   = nys - 1
+            mbeg  = jcs
+            moff  = 0
          ENDIF
       ELSEIF ( edge == 'n' )  THEN
-         j  = nyn + 1
-         jb = nyn + 2
+         IF ( var == 'v' )  THEN
+            j     = nyn + 1
+            jend  = nyn + 1
+            jawbc = 0
+            jbc   = nyn + 1
+            mbeg  = jcn - 2
+            moff  = 0
+         ELSE
+            j     = nyn + 1
+            jend  = nyn + jgsr
+            jawbc = 0
+            jbc   = nyn + 1
+            mbeg  = jcn - 2
+            moff  = 1
+         ENDIF
+      ENDIF
+
+      IF  ( var == 'w' )  THEN
+         ntopw = 1
+      ELSE
+         ntopw = 0
+      ENDIF
+
+      IF  ( var == 'u' )  THEN
+         lrightu = 0
+      ELSE
+         lrightu = 1
       ENDIF
 
@@ -5840 +7148 @@
          var_flag = 0
       ENDIF
-
-      DO  i = nxl, nxr+1
-!
-!--      Determine vertical index of topography top at grid point (j,i)
-!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)
-            n = kc(k)              
-            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
-            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
-            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
-            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
-            fk       = r1y(j) * fkj  + r2y(j) * fkjp
-            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
-            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+!AH
+!
+!--   Substitute the necessary parent-grid data to the work array workarrc_sn. 
+      workarrc_sn = 0.0_wp     
+      IF  ( pdims(1) > 1 )  THEN
+#if defined( __parallel )
+         IF  ( nxl == 0 )  THEN   ! if ( bc_dirichlet_l )
+            workarrc_sn(0:cg%nz+1,0:2,iclw:icrw-1)                              &
+                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1)
+         ELSE IF  ( nxr == nx )  THEN    ! if ( bc_dirichlet_r )
+            workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw)                              &
+                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw)
+         ELSE
+            workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                            &
+                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
+         ENDIF
+!
+!--      Left-right exchange if more than one PE subdomain in the x-direction.
+!--      Note that in case of 3-D nesting the left (pleft == MPI_PROC_NULL) and 
+!--      right (pright == MPI_PROC_NULL) boundaries are not exchanged because 
+!--      the nest domain is not cyclic.
+!--      From left to right
+         CALL MPI_SENDRECV( workarrc_sn(0,0,iclw+1), 1,                         &
+              workarrc_sn_exchange_type, pleft,   0,                            &
+              workarrc_sn(0,0,icrw), 1,                                         &
+              workarrc_sn_exchange_type, pright,  0,                            &
+              comm2d, status, ierr ) 
+!
+!--      From right to left        
+         CALL MPI_SENDRECV( workarrc_sn(0,0,icrw-1), 1,                         &
+              workarrc_sn_exchange_type, pright,  1,                            &
+              workarrc_sn(0,0,iclw), 1,                                         &
+              workarrc_sn_exchange_type, pleft,   1,                            &
+              comm2d, status, ierr )
+#endif      
+      ELSE
+         workarrc_sn(0:cg%nz+1,0:2,iclw+1:icrw-1)                               &
+                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw+1:icrw-1)
+      ENDIF
+!
+!AH
+
+      workarr_sn = 0.0_wp
+      
+      DO  i = nxl-1, nxr+1
+!AH      DO  i = nxl, nxr
+         DO  ja = j, jend
+            jaw = ja - j
+            DO  k = nzb, nzt+1
+               l   = ic(i)
+               lp1 = MIN( l + 1, icrw )  ! If l+1 > icr (l=ic(nxr+1)), r1x = 1 and r2x = 0
+!AH               m   = jc(ja) - jc(j)
+               m   = jc(ja) - mbeg
+               mp1 = MIN( m + 1, 2 )  ! If m+1 > 2 (m=jc(nyn+1)-mbeg), r1y = 1 and r2y = 0
+               n   = kc(k)
+               np1 = n + 1
+!AH
+!               fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
+!               fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
+!               fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
+!               fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
+!AH
+               fkj      = r1x(i) * workarrc_sn(n,m,l)     + r2x(i) * workarrc_sn(n,m,lp1)
+               fkjp     = r1x(i) * workarrc_sn(n,mp1,l)   + r2x(i) * workarrc_sn(n,mp1,lp1)
+               fkpj     = r1x(i) * workarrc_sn(np1,m,l)   + r2x(i) * workarrc_sn(np1,m,lp1)
+               fkpjp    = r1x(i) * workarrc_sn(np1,mp1,l) + r2x(i) * workarrc_sn(np1,mp1,lp1)
+
+               fk       = r1y(ja) * fkj  + r2y(ja) * fkjp
+               fkp      = r1y(ja) * fkpj + r2y(ja) * fkpjp
+!AH
+!               f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+               workarr_sn(k,jaw,i) = r1z(k) * fk   + r2z(k) * fkp
+
+!               if  ( ( edge == 's' ) .and. ( ja  == jbc ) ) then
+!                  write(9,"('pmci_interp_tril_sn: ',a2,2x,11(i4,2x),1(e12.5,2x))") var, k, ja, i, jbc, jaw, n, m, l, np1, mp1, lp1, &
+!                       workarr_sn(k,jaw,i)
+!                  flush(9)
+!               endif
+
+!AH               
+            ENDDO
          ENDDO
-      ENDDO
+      ENDDO      
 !
 !--   Generalized log-law-correction algorithm.
@@ -5938 +7310 @@
       ENDIF  ! ( topography /= 'flat' )
 !
-!--   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
+!--   Apply the reversibility correction. 
+      
+!      if  ( var == 'u' )  then
+
+      m  = jc(jbc) + moff
+      mw = 1
+!      write(9,"('pmci_interp_tril_sn: edge, var, m, j, jend, jfl(m), jfu(m) = ',2(a2,2x),5(i4,2x))") &
+!           edge, var, m, j, jend, jfl(m), jfu(m)
+!      flush(9)
+!AH      DO  l = iclw+1, icrw-1
+      DO  l = iclw + 1, icrw - lrightu   ! lrightu = 0 for u and 1 for all others
+         DO  n = 0, kct + ntopw            ! ntopw = 1 for w and 0 for all others
+            ijk = 1
+            cellsum   = 0.0_wp 
+            cellsumd  = 0.0_wp
+            DO  i = ifl(l), ifu(l)
+               iw = MAX( MIN( i, nx+1 ), -1 )
+               DO  ja = jfl(m), jfu(m)
+                  jaw = ja - j
+                  jw  = MAX( MIN( ja, ny+1 ), -1 )
+                  DO  k = kfl(n), kfu(n)
+                     kw = MIN( k, nzt+1 )
+!AH                     cellsum = cellsum + MERGE( f(k,ja,i), 0.0_wp,              &
+!AH                          BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+                     cellsum = cellsum + MERGE( workarr_sn(k,jaw,i), 0.0_wp,       &
+                          BTEST( wall_flags_0(kw,jw,iw), var_flag ) )                    
+!AH                     celltmpd(ijk) = ABS( fc(n,m,l) - f(k,ja,i) )
+!AH                     celltmpd(ijk) = ABS( workarrc_sn(n,mw,l) - f(k,ja,i) )
+                     celltmpd(ijk) = ABS( workarrc_sn(n,mw,l) - workarr_sn(k,jaw,i) )
+                     cellsumd      = cellsumd  + MERGE( celltmpd(ijk),          &
+                          0.0_wp, BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+
+!                     write(9,"('sn1: ',a1,2x,8(i4,2x),5(e12.5,2x))") var, n, m, l, k, ja, i, jaw, ijk,  &
+!                          workarrc_sn(n,mw,l), workarr_sn(k,jaw,i), cellsum, celltmpd(ijk), cellsumd
+!                     flush(9)
+
+                     ijk = ijk + 1
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            IF ( ijkfc(n,m,l) /= 0 )  THEN
+               cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+!AH               rcorr     = fc(n,m,l) - cellsum
+               rcorr     = workarrc_sn(n,mw,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)
+!AH                  jaw = ja - jfl(m)
+                  jaw = ja - j
+                  DO  k = kfl(n), kfu(n)
+                     rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+!AH                     f(k,ja,i) = f(k,ja,i) + rcorr_ijk
+                     workarr_sn(k,jaw,i) = workarr_sn(k,jaw,i) + rcorr_ijk
+
+!                     write(9,"('sn2: ',a1,2x,9(i4,2x),4(e12.5,2x))") var, n, m, l, k, ja, i, jaw, ijk, ijkfc(n,m,l), &
+!                          rcorr, rcorr_ijk, workarr_sn(k,jaw,i), workarrc_sn(n,mw,l)
+!                     flush(9)
+
+                     ijk = ijk + 1
+                  ENDDO
+               ENDDO
+            ENDDO
+!
+!--         Velocity components tangential to the boundary need special 
+!--         treatment because their anterpolation cells are flat in their 
+!--         direction and hence do not cover all the child-grid boundary nodes. 
+!--         These components are next linearly interpolated to those child-grid 
+!--         boundary nodes which are not covered by the anterpolation cells.     
+            IF  ( ( var == 'w' ) .AND. ( n > 0 ) )  THEN
+               DO  k = kfu(n-1)+1, kfl(n)-1
+                  IF  ( k <= nzt )  THEN
+                     DO  i = ifl(l), ifu(l)
+                        IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+!AH                         f(k,ja,i) = r1z(k) * f(kfu(n-1),ja,i)                & 
+!AH                               + r2z(k) * f(kfl(n),ja,i)
+                           workarr_sn(k,jawbc,i) =                              &
+                                r1z(k) * workarr_sn(kfu(n-1),jawbc,i) +         & 
+                                r2z(k) * workarr_sn(kfl(n),jawbc,i)
+   
+!                           write(9,"('sn3: ',a1,2x,7(i4,2x),3(e12.5,2x))") var, n, m, l, k, jawbc, i, jbc,  &
+!                                workarr_sn(k-1,jawbc,i), workarr_sn(k,jawbc,i), workarr_sn(k+1,jawbc,i)
+!                           flush(9)
+
+                        ENDIF
+                     ENDDO
+                  ENDIF
+               ENDDO
+            ENDIF
+               
+            IF  ( ( var == 'u' ) .AND. ( l > iclw ) )  THEN
+               DO  i = ifu(l-1)+1, ifl(l)-1
+                  IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+                     DO  k = kfl(n), kfu(n)
+                        IF  ( k <= nzt+1 )  THEN
+!AH                         f(k,ja,i) = r1x(i) * f(k,ja,ifu(l-1))                &
+!AH                               + r2x(i) * f(k,ja,ifl(l))
+                           workarr_sn(k,jawbc,i) =                              &
+                                r1x(i) * workarr_sn(k,jawbc,ifu(l-1)) +         &
+                                r2x(i) * workarr_sn(k,jawbc,ifl(l))
+                           
+!                           write(9,"('sn4: ',a1,2x,7(i4,2x),3(e12.5,2x))") var, n, m, l, k, jawbc, i, jbc, &
+!                                workarr_sn(k,jawbc,i-1), workarr_sn(k,jawbc,i), workarr_sn(k,jawbc,i+1)
+!                           flush(9)
+
+                        ENDIF
+                     ENDDO
+                  ENDIF
+               ENDDO
+            ENDIF
+               
+         ENDDO  ! n
+      ENDDO  ! l 
+
+!      endif  ! var
+
+!
+!--   Finally substitute the boundary values.
+      f(nzb:nzt+1,jbc,nxl:nxr) = workarr_sn(nzb:nzt+1,jawbc,nxl:nxr)
+
+!      do  k = 0, 2
+!         do  i = nxl, nxr 
+!            if  ( edge == 's' )  then
+!               write(9,"('sn5: ',2(a2,2x),4(i4,2x),4(e12.5,2x))") edge, var, k, i, jbc, jawbc,  &
+!                    workarr_sn(k,jawbc,i), f(k,jbc,i), f(k,jbc+1,i), f(k,jbc+2,i)
+!            else if  ( edge == 'n' )  then
+!               write(9,"('sn5: ',2(a2,2x),4(i4,2x),4(e12.5,2x))") edge, var, k, i, jbc, jawbc,  &
+!                    f(k,jbc-2,i), f(k,jbc-1,i), f(k,jbc,i), workarr_sn(k,jawbc,i)
+!            endif
+!         enddo
+!      enddo
+!      flush(9)
+
+   END SUBROUTINE pmci_interp_tril_sn
+
+ 
+
+   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 )
+
+!
+!--   Interpolation of ghost-node values used as the child-domain boundary
+!--   conditions. This subroutine handles the top boundary. 
+!--   This subroutine is based on trilinear interpolation.
+
+      IMPLICIT NONE
+
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
+                                      INTENT(INOUT) ::  f     !< Child-grid array
+      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr),                          &
+                                      INTENT(IN)    ::  fc    !< Parent-grid array
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r1x   !<
+      REAL(wp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  r2x   !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r1y   !<
+      REAL(wp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  r2y   !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z   !<
+!AH      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z   !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r1z   !<
+      REAL(wp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  r2z   !<
+
+      
+      INTEGER(iwp), DIMENSION(nxlfc:nxrfc), INTENT(IN)  ::  ic    !<
+      INTEGER(iwp), DIMENSION(nysfc:nynfc), INTENT(IN)  ::  jc    !<
+!AH      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  kc    !<
+      INTEGER(iwp), DIMENSION(nzb:nzt+kgsr), INTENT(IN) ::  kc    !<
+
+      INTEGER(iwp) :: kct
+!AH
+!      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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+      INTEGER(iwp), DIMENSION(jcsa:jcna), 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
+!AH
+!      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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
+
+      CHARACTER(LEN=1), INTENT(IN) :: var   !<
+
+      INTEGER(iwp) ::  i    !<
+      INTEGER(iwp) ::  ib   !<
+      INTEGER(iwp) ::  iclc !< Lower i-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  icrc !< Upper i-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  ie   !<
+      INTEGER(iwp) ::  ierr !< MPI error code
+      INTEGER(iwp) ::  ijk  !<
+      INTEGER(iwp) ::  iw   !<
+      INTEGER(iwp) ::  j    !<
+      INTEGER(iwp) ::  jb   !<
+      INTEGER(iwp) ::  jcsc !< Lower j-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  jcnc !< Upper j-index limit for copying fc-data to workarrc_t
+      INTEGER(iwp) ::  je   !< 
+      INTEGER(iwp) ::  jw   !<     
+      INTEGER(iwp) ::  k    !< Vertical child-grid index fixed to the boundary-value level 
+      INTEGER(iwp) ::  ka   !< Running vertical child-grid index
+      INTEGER(iwp) ::  kw   !< 
+      INTEGER(iwp) ::  l    !< Parent-grid index in x-direction
+      INTEGER(iwp) ::  lp1  !< l+1
+      INTEGER(iwp) ::  m    !< Parent-grid index in y-direction
+      INTEGER(iwp) ::  mp1  !< m+1
+      INTEGER(iwp) ::  n    !< Parent-grid work array index in z-direction
+      INTEGER(iwp) ::  np1  !< n+1
+      INTEGER(iwp) ::  noff !< n-offset needed on the top boundary to correctly refer to boundary ghost points  
+      INTEGER(iwp) ::  nw   !< n-index for workarrc_t
+      INTEGER(iwp) ::  var_flag  !<
+      
+      REAL(wp) ::  cellsum     !<
+      REAL(wp) ::  cellsumd    !<
+      REAL(wp) ::  fk          !<
+      REAL(wp) ::  fkj         !<
+      REAL(wp) ::  fkjp        !<
+      REAL(wp) ::  fkpj        !<
+      REAL(wp) ::  fkpjp       !<
+      REAL(wp) ::  fkp         !<
+      REAL(wp) ::  rcorr       !<
+      REAL(wp) ::  rcorr_ijk   !<
+
+      integer(iwp) :: kend
+
+
+      IF ( var == 'w' )  THEN
+         k    = nzt
+         noff = 0
+         kend = nzt
+      ELSE
+         k    = nzt + 1
+         noff = 1
+         kend = nzt+kgsr
+      ENDIF
+!
+!--   These exceedings by one are needed only to avoid stripes
+!--   and spots in visualization. They have no effect on the 
+!--   actual solution.     
+!
+!AH   These loop bounds lead to overflows with the current reduced icl, icr, jcs, jcn
+!AH      ib = nxl-1
+!AH      ie = nxr+1
+!AH      jb = nys-1
+!AH      je = nyn+1
+      ib = nxl
+      ie = nxr
+      jb = nys
+      je = nyn
+
+      IF ( var == 'u' )  THEN 
+         var_flag = 1 
+         ie       = nxr + 1 ! Needed for finishing interpolation for u 
+      ELSEIF ( var == 'v' )  THEN
+         var_flag = 2      
+         je       = nyn + 1 ! Needed for finishing interpolation for v 
+      ELSEIF ( var == 'w' )  THEN
+         var_flag = 3
+      ELSE
+         var_flag = 0
+      ENDIF
+         
+!AH      DO  i = ib, ie
+!AH         DO  j = jb, je
+!AH            l = ic(i)
+!AH            m = jc(j)
+!AH            n = kc(k)            
+!AH            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
+!AH            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
+!AH            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
+!AH            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
+!AH            fk       = r1y(j) * fkj  + r2y(j) * fkjp
+!AH            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
+!AH            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+!AH         ENDDO
+!AH      ENDDO
+
+!      write(9,"('workarrc_t kbounds: ',a2,2x,5(i3,2x))") var, k, kend, kc(k), kc(kend), cg%nz+1
+!      flush(9)
+!AH
+!
+!--   Substitute the necessary parent-grid data to the work array. 
+!--   Note that the dimension of workarrc_t is (0:2,jcsw:jcnw,iclw:icrw),
+!--   And the jc?w and ic?w-index bounds depend on the location of the PE-
+!--   subdomain relative to the side boundaries. 
+      iclc = iclw + 1
+      icrc = icrw - 1      
+      jcsc = jcsw + 1
+      jcnc = jcnw - 1
+      IF  ( bc_dirichlet_l )  THEN
+         iclc = iclw
+      ENDIF
+      IF  ( bc_dirichlet_r )  THEN
+         icrc = icrw
+      ENDIF
+      IF  ( bc_dirichlet_s )  THEN
+         jcsc = jcsw
+      ENDIF
+      IF  ( bc_dirichlet_n )  THEN
+         jcnc = jcnw
+      ENDIF
+      workarrc_t = 0.0_wp
+      workarrc_t(0:2,jcsc:jcnc,iclc:icrc)                                       &
+           = fc(kc(k):kc(k)+2,jcsc:jcnc,iclc:icrc)
+!
+!--   Left-right exchange if more than one PE subdomain in the x-direction.
+!--   Note that in case of 3-D nesting the left and right boundaries are 
+!--   not exchanged because the nest domain is not cyclic.
+#if defined( __parallel )
+      IF  ( pdims(1) > 1 )  THEN
+!
+!--      From left to right
+         CALL MPI_SENDRECV( workarrc_t(0,jcsw,iclw+1), 1,                       &
+              workarrc_t_exchange_type_y, pleft,  0,                            &
+              workarrc_t(0,jcsw,icrw), 1,                                       &
+              workarrc_t_exchange_type_y, pright, 0,                            &
+              comm2d, status, ierr ) 
+!
+!--      From right to left        
+         CALL MPI_SENDRECV( workarrc_t(0,jcsw,icrw-1), 1,                       &
+              workarrc_t_exchange_type_y, pright, 1,                            &
+              workarrc_t(0,jcsw,iclw), 1,                                       &
+              workarrc_t_exchange_type_y, pleft,  1,                            &
+              comm2d, status, ierr )
+      ENDIF
+!
+!--   South-north exchange if more than one PE subdomain in the y-direction.
+!--   Note that in case of 3-D nesting the south and north boundaries are 
+!--   not exchanged because the nest domain is not cyclic.
+      IF  ( pdims(2) > 1 )  THEN
+!
+!--      From south to north         
+         CALL MPI_SENDRECV( workarrc_t(0,jcsw+1,iclw), 1,                       &
+              workarrc_t_exchange_type_x, psouth, 2,                            &
+              workarrc_t(0,jcnw,iclw), 1,                                       &
+              workarrc_t_exchange_type_x, pnorth, 2,                            &
+              comm2d, status, ierr ) 
+!
+!--      From north to south        
+         CALL MPI_SENDRECV( workarrc_t(0,jcnw-1,iclw), 1,                       &
+              workarrc_t_exchange_type_x, pnorth, 3,                            &
+              workarrc_t(0,jcsw,iclw), 1,                                       &
+              workarrc_t_exchange_type_x, psouth, 3,                            &
+              comm2d, status, ierr )
+      ENDIF
+#endif      
+!
+!AH
+
+      workarr_t = 0.0_wp
+
+      DO  i = ib, ie
+         DO  j = jb, je
+            DO ka = k, kend
+               l    = ic(i)
+               lp1  = MIN( l + 1, icrw ) ! If l+1 > icr (l=ic(nxr+1)), r1x = 1 and r2x = 0
+               m    = jc(j)
+               mp1  = MIN( m + 1, jcnw ) ! If m+1 > jcn (m=jc(nyn+1)), r1y = 1 and r2y = 0
+!AH
+!               n    = kc(ka)            
+               n   = kc(ka) - kc(k)
+               np1 = n + 1
+!AH
+!               fkj       = r1x(i)  * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
+!               fkjp      = r1x(i)  * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
+!               fkpj      = r1x(i)  * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
+!               fkpjp     = r1x(i)  * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
+!AH
+               fkj       = r1x(i)  * workarrc_t(n,m,l)     + r2x(i) * workarrc_t(n,m,lp1)
+               fkjp      = r1x(i)  * workarrc_t(n,mp1,l)   + r2x(i) * workarrc_t(n,mp1,lp1)
+               fkpj      = r1x(i)  * workarrc_t(np1,m,l)   + r2x(i) * workarrc_t(np1,m,lp1)
+               fkpjp     = r1x(i)  * workarrc_t(np1,mp1,l) + r2x(i) * workarrc_t(np1,mp1,lp1)
+!AH
+               fk        = r1y(j)  * fkj  + r2y(j) * fkjp
+               fkp       = r1y(j)  * fkpj + r2y(j) * fkpjp
+               workarr_t(ka,j,i) = r1z(ka) * fk + r2z(ka) * fkp
+          
+            ENDDO
+         ENDDO
+      ENDDO
+
+!AH!
+!AH!--   Just fill up the redundant second ghost-node layer for w.
+!AH      IF ( var == 'w' )  THEN
+!AH         f(nzt+1,:,:) = f(nzt,:,:)
+!AH      ENDIF
+!
+!--      Apply the reversibility correction. 
+         n  = kc(k) + noff  
+         nw = 0
+!AH         DO  l = icl-1, icr+1
+         DO  l = iclw, icrw
+!AH            DO  m = jcs-1, jcn+1
+            DO  m = jcsw, jcnw
                ijk = 1
-               cellsum   = 0.0_wp 
+               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
+                  iw = MAX( MIN( i, nx+1 ), -1 )
+                  IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+                     DO  j = jfl(m), jfu(m)
+                        jw = MAX( MIN( j, ny+1 ), -1 ) 
+                        IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+                           DO  ka = kfl(n), kfu(n)
+                              kw = MIN( ka, nzt+1 )
+!AH                            cellsum = cellsum + MERGE( f(ka,j,i), 0.0_wp,           &
+                              cellsum = cellsum + MERGE( workarr_t(ka,j,i), 0.0_wp,   &
+                                   BTEST( wall_flags_0(kw,jw,iw), var_flag ) )
+!                              cellsum = cellsum + workarr_t(ka,j,i)
+!AH                            celltmpd(ijk) = ABS( fc(n,m,l) - f(ka,j,i) )
+!AH
+!                              celltmpd(ijk) = ABS( fc(n,m,l) - workarr_t(ka,j,i) )
+                              celltmpd(ijk) = ABS( workarrc_t(nw,m,l) - workarr_t(ka,j,i) )
+!AH
+                              cellsumd      = cellsumd  + MERGE( celltmpd(ijk),       &
+                                   0.0_wp, BTEST( wall_flags_0(kw,jw,iw), var_flag ) ) 
+                              ijk = ijk + 1
+                           ENDDO
+                        ENDIF
                      ENDDO
-                  ENDDO
+                  ENDIF
                ENDDO
 
                IF ( ijkfc(n,m,l) /= 0 )  THEN
                   cellsum   = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
-                  rcorr     = fc(n,m,l) - cellsum
+!AH
+!                  rcorr     = fc(n,m,l) - cellsum
+                  rcorr     = workarrc_t(nw,m,l) - cellsum
+!AH                  
                   cellsumd  = cellsumd / REAL( ijkfc(n,m,l), KIND=wp )
                ELSE
-                  cellsum   = 0.0_wp                  
+                  cellsum   = 0.0_wp
                   rcorr     = 0.0_wp
                   cellsumd  = 1.0_wp
@@ -5983 +7779 @@
                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  ! n
-         ENDDO  ! l 
-
-      ENDIF  ! var not u or w
-!
-!--   Store the boundary values also into the other redundant ghost node layers.
-!--   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
-         DO  jbgp = -nbgp, jb
-            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
-         ENDDO
-      ELSEIF ( edge == 'n' )  THEN
-         DO  jbgp = jb, ny+nbgp
-            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
-         ENDDO
-      ENDIF
-
-   END SUBROUTINE pmci_interp_tril_sn
-
- 
-
-   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 )
-
-!
-!--   Interpolation of ghost-node values used as the child-domain boundary
-!--   conditions. This subroutine handles the top boundary. 
-!--   This subroutine is based on trilinear interpolation.
-
-      IMPLICIT NONE
-
-      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), DIMENSION(nxlg:nxrg), INTENT(IN) ::  ic    !<
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) ::  jc    !<
-      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   !<
-
-      INTEGER(iwp) ::  i   !<
-      INTEGER(iwp) ::  ib  !<
-      INTEGER(iwp) ::  ie  !<
-      INTEGER(iwp) ::  ijk !<
-      INTEGER(iwp) ::  j   !<
-      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         !<
-      REAL(wp) ::  fkjp        !<
-      REAL(wp) ::  fkpj        !<
-      REAL(wp) ::  fkpjp       !<
-      REAL(wp) ::  fkp         !<
-      REAL(wp) ::  rcorr       !<
-      REAL(wp) ::  rcorr_ijk   !<
-
-
-      IF ( var == 'w' )  THEN
-         k  = nzt
-      ELSE
-         k  = nzt + 1
-      ENDIF
-!
-!--   These exceedings by one are needed only to avoid stripes
-!--   and spots in visualization. They have no effect on the 
-!--   actual solution.     
-      ib = nxl-1
-      ie = nxr+1
-      jb = nys-1
-      je = nyn+1
-
-      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 = ib, ie
-         DO  j = jb, je
-            l = ic(i)
-            m = jc(j)
-            n = kc(k)            
-            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
-            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
-            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
-            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
-            fk       = r1y(j) * fkj  + r2y(j) * fkjp
-            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
-            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
-         ENDDO
-      ENDDO
-!
-!--   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
+                  IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+                     DO  j = jfl(m), jfu(m)
+                        IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+                           DO  ka = kfl(n), kfu(n)
+                              rcorr_ijk = rcorr * celltmpd(ijk) / cellsumd 
+!AH                            f(ka,j,i) = f(ka,j,i) + rcorr_ijk
+                              workarr_t(ka,j,i) = workarr_t(ka,j,i) + rcorr_ijk
+
+!                                 if  ( i == 128 .and. var == 'u' )  then
+!                                    write(9,"('t2: ', 8(i4,2x),3(e12.5,2x))") nw, m, l, ka, j, i, ifl(l), ifu(l), &
+!                                         rcorr, rcorr_ijk, workarr_t(ka,j,i)
+!                                    flush(9)
+!                                 endif
+
+                              ijk = ijk + 1
+                           ENDDO
+                        ENDIF
                      ENDDO
+                  ENDIF
+               ENDDO
+!
+!--            Velocity components tangential to the boundary need special 
+!--            "finishing" because their anterpolation cells are flat in their 
+!--            direction and hence do not cover all the child-grid boundary nodes. 
+!--            These components are next linearly interpolated to those child-grid 
+!--            boundary nodes which are not covered by the anterpolation cells. 
+               IF  ( ( var == 'v' ) .AND. ( m > jcs ) )  THEN
+                  DO  j = jfu(m-1)+1, jfl(m)-1
+                     IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+                        DO  i = ifl(l), ifu(l)
+                           IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+                              DO  ka = kfl(n), kfu(n)        ! This loop should be removed and fixed k-value be used instead
+!AH                              f(ka,j,i) = r1y(j) * f(ka,jfu(m-1),i)             &
+!AH                                   + r2y(j) * f(ka,jfl(m),i)
+                                 workarr_t(ka,j,i) = r1y(j) * workarr_t(ka,jfu(m-1),i) &
+                                      + r2y(j) * workarr_t(ka,jfl(m),i)
+                              ENDDO
+                           ENDIF
+                        ENDDO
+                     ENDIF
                   ENDDO
-               ENDDO
+               ENDIF
+               
+               IF  ( ( var == 'u' ) .AND. ( l > icl ) )  THEN
+                  DO  i = ifu(l-1)+1, ifl(l)-1
+                     IF  ( ( i >= nxl ) .AND. ( i <= nxr+1 ) )  THEN
+                        DO  j = jfl(m), jfu(m)
+                           IF  ( ( j >= nys ) .AND. ( j <= nyn+1 ) )  THEN
+                              DO  ka = kfl(n), kfu(n)        ! This loop should be removed and fixed k-value be used instead
+!AH                              f(ka,j,i) = r1x(i) * f(ka,j,ifu(l-1)) +           &
+!AH                                   r2x(i) * f(ka,j,ifl(l))
+                                 workarr_t(ka,j,i) = r1x(i) * workarr_t(ka,j,ifu(l-1)) &
+                                      + r2x(i) * workarr_t(ka,j,ifl(l))
+
+!                                 if  ( i == 127 )  then
+!                                    write(9,"('t4: ', 8(i4,2x),3(e12.5,2x))") nw, m, l, ka, j, ifu(l-1), i, ifl(l), &
+!                                         workarr_t(ka,j,ifu(l-1)), workarr_t(ka,j,i), workarr_t(ka,j,ifl(l))
+!                                    flush(9)
+!                                 endif
+
+                              ENDDO
+                           ENDIF
+                        ENDDO
+                     ENDIF
+                  ENDDO
+               ENDIF
                
             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,:,:)
+!
+!--   Finally substitute the boundary values.
+      f(k,nys:nyn,nxl:nxr) = workarr_t(k,nys:nyn,nxl:nxr)
+      IF  ( var == 'w' )  THEN
+         f(k+1,nys:nyn,nxl:nxr) = f(k,nys:nyn,nxl:nxr)
       ENDIF
 
@@ -6204 +7874 @@
 
 !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
+!AH
+!       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), DIMENSION(0:cg%nz+1,jcsa:jcna,icla:icra), INTENT(IN) :: ijkfc !< number of child grid points contributing to a parent grid box
+!AH
 
        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) ::  iclant    !< Left boundary index for anterpolation along x
+       INTEGER(iwp) ::  icrant    !< Right boundary index for anterpolation along x
        INTEGER(iwp) ::  ii        !< Running index x-direction - coarse grid
        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) ::  jcnant    !< North boundary index for anterpolation along y
+       INTEGER(iwp) ::  jcsant    !< South boundary index for anterpolation along y
        INTEGER(iwp) ::  jj        !< Running index y-direction - coarse grid
        INTEGER(iwp) ::  k         !< Running index z-direction - fine-grid     
        INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
+       INTEGER(iwp) ::  kctant    !< Top boundary index for anterpolation along z
        INTEGER(iwp) ::  kk        !< Running index z-direction - coarse grid
        INTEGER(iwp) ::  var_flag  !< bit number used to flag topography on respective grid
@@ -6221 +7895 @@
        INTEGER(iwp), INTENT(IN) ::  kct   !< Top boundary index for anterpolation along z
 
-       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(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(icla:icra), INTENT(IN) ::  ifl !< Indicates start index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(icla:icra), INTENT(IN) ::  ifu !< Indicates end index of child cells belonging to certain parent cell - x direction
+       INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfl !< Indicates start index of child cells belonging to certain parent cell - y direction
+       INTEGER(iwp), DIMENSION(jcsa:jcna), INTENT(IN) ::  jfu !< Indicates start index of child cells belonging to certain parent cell - y direction
+!AH
+
 !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
@@ -6237 +7918 @@
   
 !
-!--    Define the index bounds icla, icra, jcsa and jcna.
+!--    Define the index bounds iclant, icrant, jcsant and jcnant.
 !--    Note that kcb is simply zero and kct enters here as a parameter and it is
 !--    determined in pmci_init_anterp_tophat.
 !--    Please note, grid points used also for interpolation (from parent to
 !--    child) are excluded in anterpolation, e.g. anterpolation is only from 
-!--    nzb:kct-1, as kct is used for interpolation. Following this approach
-!--    avoids numerical problems which may accumulate, particularly for shallow
-!--    child domain, leading to increased velocity variances. A more 
-!--    comprehensive explanation for this is still pending.
-       icla = coarse_bound_anterp(1)
-       icra = coarse_bound_anterp(2) 
-       jcsa = coarse_bound_anterp(3)
-       jcna = coarse_bound_anterp(4)
+!--    nzb:kct-1, as kct is used for interpolation.
+       iclant = icl
+       icrant = icr
+       jcsant = jcs
+       jcnant = jcn
+       kctant = kct - 1
+
        kcb  = 0
        IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
-             IF ( var == 'u' )  THEN
-                icla = coarse_bound_anterp(1) + 2
-             ELSE
-                icla = coarse_bound_anterp(1) + 1
-             ENDIF
+             iclant = icl + 3
           ENDIF
           IF ( bc_dirichlet_r )  THEN
-             icra = coarse_bound_anterp(2) - 1
+!             icrant = icr - 5
+             icrant = icr - 3
           ENDIF
 
           IF ( bc_dirichlet_s )  THEN
-             IF ( var == 'v' )  THEN
-                jcsa = coarse_bound_anterp(3) + 2
-             ELSE
-                jcsa = coarse_bound_anterp(3) + 1
-             ENDIF
+             jcsant = jcs + 3
           ENDIF
           IF ( bc_dirichlet_n )  THEN
-             jcna = coarse_bound_anterp(4) - 1
+             jcnant = jcn - 3
           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 
@@ -6292 +7961 @@
 !--    Note that ii, jj, and kk are coarse-grid indices and i,j, and k 
 !--    are fine-grid indices.
-       DO  ii = icla, icra
-          DO  jj = jcsa, jcna
+       DO  ii = iclant, icrant
+          DO  jj = jcsant, jcnant
 !
 !--          For simplicity anterpolate within buildings and under elevated
 !--          terrain too
-             DO  kk = kcb, kct - 1                
+             DO  kk = kcb, kctant !kct - 1
                 cellsum = 0.0_wp
                 DO  i = ifl(ii), ifu(ii)
@@ -6310 +7979 @@
 !--             Spatial under-relaxation.
 !--             The relaxation buffer zones are no longer needed with
-!--             the new reversible interpolation algorithm. 23.19.2018.
+!--             the new reversible interpolation algorithm. 23.10.2018.
 !                fra  = frax(ii) * fray(jj) * fraz(kk)                
 !
@@ -6318 +7987 @@
 !--             particular for the temperature. Therefore, in case cellsum is 
 !--             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 )
+!AH                   fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +         &
+!AH                        fra * cellsum  /                                    &
+!AH                        REAL( ijkfc(kk,jj,ii), KIND=wp )
                    fc(kk,jj,ii) = cellsum / REAL( ijkfc(kk,jj,ii), KIND=wp )
                 ENDIF
 
              ENDDO
-
           ENDDO
        ENDDO

palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header

@@ -37 +37 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_cbm4/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstatp/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsagas/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90

@@ -41 +41 @@
 ! PALM. If not,see <http://www.gnu.org/licenses/>.
 !
-! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 1997-2019 Leibniz Universitaet Hannover
 !--------------------------------------------------------------------------------!
 !

palm/trunk/UTIL/inifor/tests/test-input-files.f90

@@ -21 +21 @@
 ! Current revisions:
 ! -----------------
-! Switched to get_datetime_file_list()
 ! 
 ! 
@@ -27 +26 @@
 ! -----------------
 ! $Id$
+! Switched to get_datetime_file_list()
+! 
+! 
+! 3678 2019-01-17 14:12:17Z eckhard
 ! Prefixed all INIFOR modules with inifor_
 !