Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3706)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3708)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Checks for parent / child grid line matching introduced.
+! Interpolation of nest-boundary-tangential velocity components revised.
+! 
+! 3697 2019-01-24 17:16:13Z hellstea
 ! Bugfix: upper k-bound in the child initialization interpolation 
 ! pmci_interp_1sto_all corrected.
@@ -1475,4 +1479,7 @@
 !--    included in the interpolation algorithms.
        CALL pmci_init_anterp_tophat
+!
+!--    Check that the child and parent grid lines match 
+       CALL pmci_check_grid_matching 
 
     ENDIF
@@ -1928,5 +1935,5 @@
             CEILING( parentdzmax / dzmin )
        ALLOCATE( celltmpd(1:acsize) )
-      
+
     END SUBROUTINE pmci_init_interp_tril
     
@@ -1941,5 +1948,5 @@
        igsr = NINT( cg%dx / dx, iwp )
        jgsr = NINT( cg%dy / dy, iwp )
-       kgsr = NINT( cg%dzu(cg%nz+1) / dzu(nzt+1), iwp )
+       kgsr = NINT( cg%dzw(1) / dzw(1), iwp )
        write(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
        flush(9)
@@ -3899,5 +3906,77 @@
     END SUBROUTINE pmci_init_anterp_tophat
 
-#endif
+
+
+    SUBROUTINE pmci_check_grid_matching 
+!
+!--    Check that the grid lines of child and parent do match  
+       IMPLICIT NONE
+       REAL(wp), PARAMETER :: tolefac = 1.0E-6_wp
+       REAL(wp) :: tolex
+       REAL(wp) :: toley
+       REAL(wp) :: tolez
+       INTEGER(iwp) :: non_matching_lower_left_corner
+       INTEGER(iwp) :: non_int_gsr_x
+       INTEGER(iwp) :: non_int_gsr_y
+       INTEGER(iwp) :: non_int_gsr_z
+
+
+       non_matching_lower_left_corner = 0
+       non_int_gsr_x = 0
+       non_int_gsr_y = 0
+       non_int_gsr_z = 0
+
+       IF  ( myid == 0 )  THEN
+
+          tolex = tolefac * dx
+          toley = tolefac * dy
+          tolez = tolefac * MINVAL( dzw )
+!
+!--       First check that the child grid lower left corner matches the paren grid lines.
+          IF  ( MOD( lower_left_coord_x, cg%dx ) > tolex ) non_matching_lower_left_corner = 1
+          IF  ( MOD( lower_left_coord_y, cg%dy ) > toley ) non_matching_lower_left_corner = 1
+!
+!--       Then check that the grid-spacing ratios in each direction are integer valued.    
+          IF  ( MOD( cg%dx, dx ) > tolex )  non_int_gsr_x = 1
+          IF  ( MOD( cg%dy, dy ) > toley )  non_int_gsr_y = 1
+          DO  n = 0, kctw+1
+             IF  ( ABS( cg%zw(n) - zw(kflw(n)) ) > tolez )  non_int_gsr_z = 1
+          ENDDO
+
+       ENDIF
+
+       CALL MPI_BCAST( non_matching_lower_left_corner, 1, MPI_INTEGER, 0,       &
+            comm2d, ierr )
+       IF  ( non_matching_lower_left_corner > 0 )  THEN
+          WRITE ( message_string, * )  'nested child domain lower left ',       &
+               'corner must match its parent grid lines'
+          CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
+       ENDIF
+
+       CALL MPI_BCAST( non_int_gsr_x, 1, MPI_INTEGER, 0, comm2d, ierr )
+       IF  ( non_int_gsr_x > 0 )  THEN
+          WRITE ( message_string, * )  'nesting grid-spacing ratio ',           &
+               '( parent dx / child dx ) must have an integer value'
+          CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
+       ENDIF
+
+       CALL MPI_BCAST( non_int_gsr_y, 1, MPI_INTEGER, 0, comm2d, ierr )
+       IF  ( non_int_gsr_y > 0 )  THEN
+          WRITE ( message_string, * )  'nesting grid-spacing ratio ',           &
+               '( parent dy / child dy ) must have an integer value'
+          CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
+       ENDIF
+
+       CALL MPI_BCAST( non_int_gsr_z, 1, MPI_INTEGER, 0, comm2d, ierr )
+       IF  ( non_int_gsr_z > 0 )  THEN
+          WRITE ( message_string, * )  'nesting grid-spacing ratio ',           &
+               '( parent dz / child dz ) must have an integer value for each z-level'
+          CALL message( 'pmci_check_grid_matching', 'PA0416', 3, 2, 0, 6, 0 )
+       ENDIF
+
+     END SUBROUTINE pmci_check_grid_matching
+
+
+#endif  
  END SUBROUTINE pmci_setup_child
 
@@ -4557,5 +4636,5 @@
              DO  m = mb, me
                 DO n = 0, kct + 1  
-               
+
                    DO i = ifl(l), ifu(l)
                       DO  j = jfl(m), jfl(m+1)-1
@@ -4594,7 +4673,7 @@
              DO  m = mb, me
                 DO n = 0, kct + 1
-                   
+                                      
                    DO i = ifl(l), ifu(l)
-                      DO  j = jfl(m), jfu(m)
+                      DO  j = jfl(m), jfu(m)                         
                          DO  k = kfl(n), MIN( kfu(n), nzt+1 )
                             f(k,j,i) = fc(n,m,l)
@@ -5785,4 +5864,8 @@
       REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
       REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+
+      real(wp), parameter :: c1 =  2.0_wp / 6.0_wp
+      real(wp), parameter :: c2 =  5.0_wp / 6.0_wp
+      real(wp), parameter :: c3 = -1.0_wp / 6.0_wp
 
 ! 
@@ -5858,86 +5941,93 @@
       ENDIF
 
-      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)
+      IF  ( var == 'v' .OR. var == 'w' )  THEN
+!
+!--      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
+
+         IF  ( var == 'v' )  THEN
+         
+            DO  m = jcsw+1, jcnw
+               DO n = 0, kct              
+!
+!--               First substitute only the matching-node values                                     
+                  f(kfl(n):kfu(n),jfl(m),ibc) = workarrc_lr(n,m,lw) 
+                  
+               ENDDO
+            ENDDO
+
+            DO  m = jcsw+1, jcnw-1
+               DO n = 0, kct 
+!
+!--               Then fill up the nodes in between with the averages
+                  DO  j = jfl(m)+1, jfl(m+1)-1
+                     f(kfl(n):kfu(n),j,ibc) = 0.5_wp * ( f(kfl(n):kfu(n),jfl(m),ibc) &
+                          + f(kfl(n):kfu(n),jfl(m+1),ibc) )
                   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) 
+
+         ELSE IF  ( var == 'w' )  THEN
+
+            DO  m = jcsw+1, jcnw-1
+               DO n = 0, kct + 1   ! It is important to go up to kct+1  
+!
+!--               First substitute only the matching-node values
+                  f(kfu(n),jfl(m):jfu(m),ibc) = workarrc_lr(n,m,lw)
+                  
+               ENDDO
+            ENDDO
+
+            DO  m = jcsw+1, jcnw-1
+               DO n = 1, kct + 1   ! It is important to go up to kct+1  
+!
+!--               Then fill up the nodes in between with the averages                  
+                  DO  k = kfu(n-1)+1, kfu(n)-1 
+                     f(k,jfl(m):jfu(m),ibc) = 0.5_wp * ( f(kfu(n-1),jfl(m):jfu(m),ibc) &
+                          + f(kfu(n),jfl(m):jfu(m),ibc) )
                   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
+
+         ENDIF
+
+      ELSE   ! u or scalars
          
          DO  m = jcsw+1, jcnw-1
@@ -6081,4 +6171,8 @@
       REAL(wp) ::  rcorr        !< Average reversibility correction for the whole anterpolation cell
       REAL(wp) ::  rcorr_ijk    !< Reversibility correction distributed to the individual child-grid nodes
+
+      real(wp), parameter :: c1 =  2.0_wp / 6.0_wp
+      real(wp), parameter :: c2 =  5.0_wp / 6.0_wp
+      real(wp), parameter :: c3 = -1.0_wp / 6.0_wp
 
 !
@@ -6154,90 +6248,93 @@
       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)
+      IF  ( var == 'u' .OR. var == 'w' )  THEN
+!
+!--      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
+
+         IF  ( var == 'u' )  THEN
+         
+            DO  l = iclw+1, icrw
+               DO n = 0, kct 
+!
+!--               First substitute only the matching-node values                   
+                  f(kfl(n):kfu(n),jbc,ifl(l)) = workarrc_sn(n,mw,l)
+
+               ENDDO
+            ENDDO
+
+            DO  l = iclw+1, icrw-1
+               DO n = 0, kct 
+!
+!--               Then fill up the nodes in between with the averages
+                  DO  i = ifl(l)+1, ifl(l+1)-1
+                     f(kfl(n):kfu(n),jbc,i) = 0.5_wp * ( f(kfl(n):kfu(n),jbc,ifl(l)) &
+                          + f(kfl(n):kfu(n),jbc,ifl(l+1)) )
                   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) 
-!AH
-!                     write(9,"('sn u: ',6(i3,2x),2(e12.5,2x))") n, m, l, k, jbc, i, fc(n,m,l), f(k,jbc,i)
-!                     flush(9)
-!AH
+
+         ELSE IF  ( var == 'w' )  THEN
+
+            DO  l = iclw+1, icrw-1
+               DO n = 0, kct + 1   ! It is important to go up to kct+1  
+!                  
+!--               First substitute only the matching-node values                 
+                  f(kfu(n),jbc,ifl(l):ifu(l)) = workarrc_sn(n,mw,l)
+
+               ENDDO
+            ENDDO
+
+            DO  l = iclw+1, icrw-1
+               DO n = 1, kct + 1   ! It is important to go up to kct+1  
+!
+!--               Then fill up the nodes in between with the averages
+                  DO  k = kfu(n-1)+1, kfu(n)-1 
+                     f(k,jbc,ifl(l):ifu(l)) = 0.5_wp * ( f(kfu(n-1),jbc,ifl(l):ifu(l)) &
+                          + f(kfu(n),jbc,ifl(l):ifu(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
+
+         ENDIF
+
+      ELSE   ! v or scalars
          
          DO  l = iclw+1, icrw-1
@@ -6358,6 +6455,11 @@
       REAL(wp) ::  fkpjp       !<
       REAL(wp) ::  fkp         !<
+      REAL(wp) ::  f_interp    !<
       REAL(wp) ::  rcorr       !<
       REAL(wp) ::  rcorr_ijk   !<
+
+      real(wp), parameter :: c1 =  2.0_wp / 6.0_wp
+      real(wp), parameter :: c2 =  5.0_wp / 6.0_wp
+      real(wp), parameter :: c3 = -1.0_wp / 6.0_wp
 
 
@@ -6380,81 +6482,192 @@
       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
-
+      nw = noff 
+!      write(9,"(a1,2x,5(i3,2x))") var, k, kc(k), noff, n, nw 
+!      flush(9)
+!
+!--   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      
+
+      IF  (  var == 'w' )  THEN 
+!         DO  l = iclw+1, icrw-1
+!            DO  m = jcsw+1, jcnw-1
+         DO  l = iclw, icrw   !  Also fill up beyond edges
+            DO  m = jcsw, jcnw   !  Also fill up beyond edges
+
+               DO  i = ifl(l), ifu(l)
+                  DO  j = jfl(m), jfu(m)
+                     f(k,j,i) = workarrc_t(nw,m,l)
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+      ENDIF
+
+      IF  ( var == 'u' )  THEN
+         
+         DO  l = iclw+1, icrw  ! Note that we have to go up to icrw here. 
+            DO  m = jcsw+1, jcnw-1
+!
+!--            1st-order upwind interpolation
+              f_interp = workarrc_t(nw,m,l)
+!
+!--            3rd-order upwind biased interpolation. Note that w is available 
+!--            here since it is already interpolated.
+!               IF  ( w(k,jfl(m),ifl(l)) + w(k,jfl(m),ifl(l)-1) < 0.0_wp )  THEN
+!                  f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l)    &
+!                       + c3 * workarrc_t(2,m,l) 
+!               ELSE
+!                  f_interp = workarrc_t(nw,m,l)
+!               ENDIF
+
+!               DO  i = ifl(l), ifl(l+1)-1
+!                  DO  j = jfl(m), jfu(m) 
+!                     f(k,j,i) = fc(n,m,l)
+!                  ENDDO
+!               ENDDO
+!
+!--            First substitute only the matching-node values 
+               f(k,jfl(m):jfu(m),ifl(l)) = f_interp
+
+            ENDDO
+         ENDDO
+!
+!--      Then fill up the nodes in between with the averages
          DO  l = iclw+1, icrw-1
             DO  m = jcsw+1, jcnw-1
-
+                  
+               DO  i = ifl(l)+1, ifl(l+1)-1
+                  f(k,jfl(m):jfu(m),i) = 0.5_wp * ( f(k,jfl(m):jfu(m),ifl(l))   &
+                       + f(k,jfl(m):jfu(m),ifl(l+1)) )
+               ENDDO
+                  
+            ENDDO
+         ENDDO
+
+      ENDIF
+      
+      IF  ( var == 'v' )  THEN
+
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw  ! Note that we have to go up to jcnw here. 
+!
+!--            1st-order upwind interpolation
+               f_interp  = workarrc_t(nw,m,l)
+!
+!--            3rd-order upwind biased interpolation. Note that w is available 
+!--            here since it is already interpolated.
+!               IF  ( w(k,jfl(m),ifl(l)) + w(k,jfl(m)-1,ifl(l)) < 0.0_wp )  THEN
+!                  f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) + c3 * workarrc_t(2,m,l) 
+!               ELSE
+!                  f_interp = workarrc_t(nw,m,l)
+!               ENDIF
+
+!               DO  i = ifl(l), ifu(l)
+!                  DO  j = jfl(m), jfl(m+1)-1
+!                     f(k,j,i) = fc(n,m,l)
+!                  ENDDO
+!               ENDDO
+!
+!--            First substitute only the matching-node values 
+               f(k,jfl(m),ifl(l):ifu(l)) = f_interp
+
+            ENDDO
+         ENDDO
+!
+!--      Then fill up the nodes in between with the averages
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+                  
+               DO  j = jfl(m)+1, jfl(m+1)-1
+                  f(k,j,ifl(l):ifu(l)) = 0.5_wp * ( f(k,jfl(m),ifl(l):ifu(l))   &
+                       + f(k,jfl(m+1),ifl(l):ifu(l)) )
+               ENDDO
+               
+            ENDDO
+         ENDDO
+
+      ENDIF
+
+      IF  (  var == 's' )  THEN 
+         DO  l = iclw+1, icrw-1
+            DO  m = jcsw+1, jcnw-1
+!
+!--            1st-order upwind interpolation
+               f_interp  = workarrc_t(nw,m,l)
+!
+!--            3rd-order upwind biased interpolation. Note that w is available 
+!--            here since it is already interpolated.
+!               IF  ( w(k,jfl(m),ifl(l)) < 0.0_wp )  THEN
+!                  f_interp = c1 * workarrc_t(0,m,l) + c2 * workarrc_t(1,m,l) + c3 * workarrc_t(2,m,l) 
+!               ELSE
+!                  f_interp = workarrc_t(nw,m,l)
+!               ENDIF
+              
                DO  i = ifl(l), ifu(l)
                   DO  j = jfl(m), jfu(m)
-                     f(k,j,i) = fc(n,m,l)
+                     f(k,j,i) = f_interp
                   ENDDO
                ENDDO
@@ -6462,48 +6675,5 @@
             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
+      ENDIF
 !
 !--   Just fill up the redundant second ghost-node layer in case of var == w.
