Index: /palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3740)
+++ /palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3741)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Interpolations and child initialization adjusted to handle set ups with child 
+! pe-subdomain dimension not integer divisible by the grid-spacing ratio in the 
+! respective direction. Set ups with pe-subdomain dimension smaller than the 
+! grid-spacing ratio in the respective direction are now forbidden.
+! 
+! 3708 2019-01-30 12:58:13Z hellstea
 ! Checks for parent / child grid line matching introduced.
 ! Interpolation of nest-boundary-tangential velocity components revised.
@@ -1504,8 +1510,4 @@
        INTEGER(iwp) :: jauxs    !<
        INTEGER(iwp) :: jauxn    !< 
-       REAL(wp) ::  loffset     !<
-       REAL(wp) ::  noffset     !<
-       REAL(wp) ::  roffset     !<
-       REAL(wp) ::  soffset     !<
        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 
@@ -1580,6 +1582,5 @@
 !--    Left
        IF  ( bc_dirichlet_l )  THEN
-          loffset = MOD( coord_x(nxl), cg%dx )
-          xexl  = 2 * cg%dx + loffset
+          xexl  = 2 * cg%dx
           iauxl = 0
        ELSE
@@ -1597,6 +1598,5 @@
 !--    Right
        IF  ( bc_dirichlet_r )  THEN
-          roffset = MOD( coord_x(nxr+1), cg%dx )
-          xexr  = 2 * cg%dx + roffset
+          xexr  = 2 * cg%dx
           iauxr = 0  
        ELSE
@@ -1607,5 +1607,5 @@
        DO  i = cg%nx, 0 , -1
           IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= xce )  THEN
-             icr = MIN( cg%nx, i )
+             icr = MIN( cg%nx, MAX( icl, i ) )
              EXIT
           ENDIF
@@ -1614,6 +1614,5 @@
 !--    South
        IF  ( bc_dirichlet_s )  THEN
-          soffset = MOD( coord_y(nys), cg%dy )
-          yexs  = 2 * cg%dy + soffset
+          yexs  = 2 * cg%dy
           jauxs = 0  
        ELSE
@@ -1631,6 +1630,5 @@
 !--    North
        IF  ( bc_dirichlet_n )  THEN
-          noffset = MOD( coord_y(nyn+1), cg%dy )
-          yexn  = 2 * cg%dy + noffset
+          yexn  = 2 * cg%dy
           jauxn = 0
        ELSE
@@ -1641,5 +1639,5 @@
        DO  j = cg%ny, 0 , -1
           IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= yce )  THEN
-             jcn = MIN( cg%ny, j )
+             jcn = MIN( cg%ny, MAX( jcs, j ) )
              EXIT
           ENDIF
@@ -3556,43 +3554,16 @@
 !
 !--    i-indices of u for each ii-index value
-!--    ii=icr is redundant for anterpolation
-       tolerance = 0.000001_wp * dx
        istart = nxlg
-!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 
-!--       use only the local k,j-child-grid plane for the anterpolation, 
-!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
+       DO  ii = icla, icra
+!
+!--       The parent and child grid lines do always match in x, hence we 
+!--       use only the local k,j-child-grid plane for the anterpolation.
           i = istart
-          DO WHILE ( coord_x(i) < cg%coord_x(ii) - tolerance  .AND.  i < nxrg )
+          DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg )
              i = i + 1
           ENDDO
-          IF ( ABS( coord_x(i) - cg%coord_x(ii) ) < tolerance )  THEN
-             i = istart
-             DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg )
-                i = i + 1
-             ENDDO
-             iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
-             ifuu(ii) = iflu(ii)
-             istart   = iflu(ii)
-          ELSE
-             i = istart
-             DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )       &
-                  .AND.  ( i < nxrg ) )
-                i  = i + 1
-             ENDDO
-             iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
-             ir = i
-             DO  WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx )     &
-                  .AND.  ( i < nxrg+1 ) )
-                i  = i + 1
-                ir = MIN( i, nxrg )
-             ENDDO
-             ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg )
-             istart = iflu(ii)
-          ENDIF
+          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
+          ifuu(ii) = iflu(ii)
+          istart   = iflu(ii)
 !
 !--       Print out the index bounds for checking and debugging purposes
@@ -3605,8 +3576,5 @@
 !
 !--    i-indices of others for each ii-index value
-!--    ii=icr is redundant for anterpolation
        istart = nxlg
-!AH       DO  ii = icl, icr
-!AH       DO  ii = icl-1, icr+1
        DO  ii = icla, icra
           i = istart
@@ -3633,42 +3601,16 @@
 !
 !--    j-indices of v for each jj-index value
-!--    jj=jcn is redundant for anterpolation
-       tolerance = 0.000001_wp * dy
        jstart = nysg
-!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 
-!--       use only the local k,i-child-grid plane for the anterpolation, 
-!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
+!--       The parent and child grid lines do always match in y, hence we 
+!--       use only the local k,i-child-grid plane for the anterpolation.
           j = jstart
-          DO WHILE ( coord_y(j) < cg%coord_y(jj) - tolerance  .AND.  j < nyng )
+          DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng )
              j = j + 1
           ENDDO
-          IF ( ABS( coord_y(j) - cg%coord_y(jj) ) < tolerance )  THEN
-             j = jstart
-             DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng )
-                j = j + 1
-             ENDDO
-             jflv(jj) = MIN( MAX( j, nysg ), nyng )
-             jfuv(jj) = jflv(jj)
-             jstart   = jflv(jj)
-          ELSE
-             j = jstart
-             DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )       &
-                  .AND.  ( j < nyng ) )
-                j  = j + 1
-             ENDDO
-             jflv(jj) = MIN( MAX( j, nysg ), nyng )
-             jr = j
-             DO  WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy )     &
-                  .AND.  ( j < nyng+1 ) )
-                j  = j + 1
-                jr = MIN( j, nyng )
-             ENDDO
-             jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng )
-             jstart = jflv(jj)
-          ENDIF
+          jflv(jj) = MIN( MAX( j, nysg ), nyng )
+          jfuv(jj) = jflv(jj)
+          jstart   = jflv(jj)
 !
 !--       Print out the index bounds for checking and debugging purposes
@@ -3680,8 +3622,5 @@
 !
 !--    j-indices of others for each jj-index value
-!--    jj=jcn is redundant for anterpolation
        jstart = nysg
-!AH       DO  jj = jcs, jcn
-!AH       DO  jj = jcs-1, jcn+1
        DO  jj = jcsa, jcna
           j = jstart
@@ -3709,44 +3648,21 @@
 !--    k-indices of w for each kk-index value
 !--    Note that anterpolation index limits are needed also for the top boundary 
-!--    ghost cell level because of the reversibility correction in the interpolation.
+!--    ghost cell level because they are used also in the interpolation.
        kstart  = 0
        kflw(0) = 0
        kfuw(0) = 0
-       tolerance = 0.000001_wp * dzw(1)
        DO kk = 1, cg%nz+1
 !
-!--       In case the child and parent grid lines match in z
-!--       use only the local j,i-child-grid plane for the anterpolation, 
-!--       i.e use 2-D anterpolation. Else, use a 3-D anterpolation.
+!--       The parent and child grid lines do always match in z, hence we 
+!--       use only the local j,i-child-grid plane for the anterpolation.
           k = kstart
-          DO WHILE ( zw(k) < cg%zw(kk) - tolerance  .AND.  k < nzt+1 )
+          DO WHILE ( ( zw(k) < cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
              k = k + 1
           ENDDO
-          IF ( ABS( zw(k) - cg%zw(kk) ) < tolerance )  THEN
-             k = kstart
-             DO WHILE ( ( zw(k) < cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
-                k = k + 1
-             ENDDO
-             kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-             kfuw(kk) = kflw(kk)
-             kstart   = kflw(kk)
-          ELSE
-             k = kstart
-             DO  WHILE ( ( zw(k) < cg%zu(kk) )  .AND.  ( k < nzt ) )
-                k = k + 1
-             ENDDO
-             kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-             IF ( kk+1 <= cg%nz+1 )  THEN
-                DO  WHILE ( ( zw(k) <= cg%zu(kk+1) )  .AND.  ( k < nzt+1 ) )
-                   k  = k + 1
-                   IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
-                ENDDO
-                kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 )
-             ELSE
-                kfuw(kk) = kflw(kk)
-             ENDIF
-             kstart = kflw(kk)
-          ENDIF
-!AH
+          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+          kfuw(kk) = kflw(kk)
+          kstart   = kflw(kk)
+!
+!--       Print out the index bounds for checking and debugging purposes
           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)
@@ -3761,5 +3677,5 @@
 !
 !--    Note that anterpolation index limits are needed also for the top boundary 
-!--    ghost cell level because of the reversibility correction in the interpolation.
+!--    ghost cell level because they are used also in the interpolation.
        DO  kk = 1, cg%nz+1
           k = kstart
@@ -3782,5 +3698,6 @@
        kflo(cg%nz+1) = nzt+kgsr ! Obsolete
        kfuo(cg%nz+1) = nzt+kgsr ! Obsolete 
-
+!
+!--       Print out the index bounds for checking and debugging purposes
        DO  kk = 1, cg%nz+1
           write(9,"('pmci_init_anterp_tophat, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))") &
@@ -3789,6 +3706,4 @@
        ENDDO
        write(9,*)
-!AH
-
 !
 !--    Precomputation of number of fine-grid nodes inside parent-grid cells.
@@ -3910,14 +3825,20 @@
     SUBROUTINE pmci_check_grid_matching 
 !
-!--    Check that the grid lines of child and parent do match  
+!--    Check that the grid lines of child and parent do match.
+!--    Also check that the child subdomain width is not smaller than 
+!--    the parent grid spacing in the respective direction. 
        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
+       REAL(wp), PARAMETER :: tolefac = 1.0E-6_wp      !< Relative tolerence for grid-line matching
+       REAL(wp) :: pesdwx                              !< Child subdomain width in x-direction 
+       REAL(wp) :: pesdwy                              !< Child subdomain width in y-direction 
+       REAL(wp) :: tolex                               !< Tolerance for grid-line matching in x-direction
+       REAL(wp) :: toley                               !< Tolerance for grid-line matching in y-direction
+       REAL(wp) :: tolez                               !< Tolerance for grid-line matching in z-direction
+       INTEGER(iwp) :: non_matching_lower_left_corner  !< Flag for non-matching lower left corner
+       INTEGER(iwp) :: non_int_gsr_x                   !< Flag for non-integer grid-spacing ration in x-direction
+       INTEGER(iwp) :: non_int_gsr_y                   !< Flag for non-integer grid-spacing ration in y-direction
+       INTEGER(iwp) :: non_int_gsr_z                   !< Flag for non-integer grid-spacing ration in z-direction
+       INTEGER(iwp) :: too_narrow_pesd_x               !< Flag for too narrow pe-subdomain in x-direction
+       INTEGER(iwp) :: too_narrow_pesd_y               !< Flag for too narrow pe-subdomain in y-direction
 
 
@@ -3926,4 +3847,6 @@
        non_int_gsr_y = 0
        non_int_gsr_z = 0
+       too_narrow_pesd_x = 0
+       too_narrow_pesd_y = 0
 
        IF  ( myid == 0 )  THEN
@@ -3944,4 +3867,9 @@
           ENDDO
 
+          pesdwx = REAL( nxr - nxl + 1, KIND=wp )
+          IF  ( pesdwx / REAL( igsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_x = 1
+          pesdwy = REAL( nyn - nys + 1, KIND=wp )
+          IF  ( pesdwy / REAL( jgsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_y = 1
+                          
        ENDIF
 
@@ -3973,5 +3901,21 @@
                '( 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
+       ENDIF    
+
+       CALL MPI_BCAST( too_narrow_pesd_x, 1, MPI_INTEGER, 0, comm2d, ierr )
+       IF  ( too_narrow_pesd_x > 0 )  THEN
+          WRITE ( message_string, * )  'child subdomain width in x-direction ',    &
+               'must not be smaller than its parent grid dx. Change the PE-grid ', &
+               'setting (npex, npey) to satisfy this requirement.'  
+          CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
+       ENDIF
+ 
+       CALL MPI_BCAST( too_narrow_pesd_y, 1, MPI_INTEGER, 0, comm2d, ierr )
+       IF  ( too_narrow_pesd_y > 0 )  THEN
+          WRITE ( message_string, * )  'child subdomain width in y-direction ',    &
+               'must not be smaller than its parent grid dy. Change the PE-grid ', &
+               'setting (npex, npey) to satisfy this requirement.'  
+          CALL message( 'pmci_check_grid_matching', 'PA0587', 3, 2, 0, 6, 0 )
+       ENDIF       
 
      END SUBROUTINE pmci_check_grid_matching
@@ -4521,6 +4465,4 @@
 !
 !--    Interpolation of the internal values for the child-domain initialization
-!--    This subroutine is based on trilinear interpolation.
-!--    Coding based on interp_tril_lr/sn/t
        IMPLICIT NONE
 
@@ -4558,5 +4500,13 @@
 
        INTEGER(iwp) ::  i        !<
+       INTEGER(iwp) ::  ib       !<
+       INTEGER(iwp) ::  ie       !<
+       INTEGER(iwp) ::  ifirst   !<
+       INTEGER(iwp) ::  ilast    !<
        INTEGER(iwp) ::  j        !<
+       INTEGER(iwp) ::  jb       !<
+       INTEGER(iwp) ::  je       !<
+       INTEGER(iwp) ::  jfirst   !<
+       INTEGER(iwp) ::  jlast    !<
        INTEGER(iwp) ::  k        !<
        INTEGER(iwp) ::  l        !<
@@ -4574,27 +4524,37 @@
        mb = jcs
        me = jcn
+       ifirst = nxl
+       ilast  = nxr
+       jfirst = nys
+       jlast  = nyn
 
        IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
-             lb = icl + 1
+             lb     = icl + 1
+             ifirst = nxl - 1
 !
 !--          For u, nxl is a ghost node, but not for the other variables
              IF ( var == 'u' )  THEN
-                lb = icl + 2
+                lb     = icl + 2
+                ifirst = nxl
              ENDIF
           ENDIF
           IF ( bc_dirichlet_s )  THEN
-             mb = jcs + 1
+             mb     = jcs + 1
+             jfirst = nys - 1 
 !
 !--          For v, nys is a ghost node, but not for the other variables
              IF ( var == 'v' )  THEN
-                mb = jcs + 2
+                mb     = jcs + 2
+                jfirst = nys
              ENDIF
           ENDIF
           IF ( bc_dirichlet_r )  THEN
-             le = icr - 1
+             le    = icr - 1
+             ilast = nxr + 1 
           ENDIF
           IF ( bc_dirichlet_n )  THEN
-             me = jcn - 1
+             me    = jcn - 1
+             jlast = nyn + 1
           ENDIF
        ENDIF
@@ -4615,4 +4575,8 @@
        IF  ( var == 'u' )  THEN
 
+          ib = ifl(lb) 
+          ie = ifl(le+1) - 1
+          jb = jfl(mb)
+          je = jfu(me)
           DO  l = lb, le
              DO  m = mb, me
@@ -4633,4 +4597,8 @@
        ELSE IF  ( var == 'v' )  THEN
 
+          ib = ifl(lb) 
+          ie = ifu(le)
+          jb = jfl(mb)
+          je = jfl(me+1) - 1
           DO  l = lb, le
              DO  m = mb, me
@@ -4651,4 +4619,8 @@
        ELSE IF  ( var == 'w' )  THEN
 
+          ib = ifl(lb) 
+          ie = ifu(le)
+          jb = jfl(mb)
+          je = jfu(me)
           DO  l = lb, le
              DO  m = mb, me
@@ -4670,4 +4642,8 @@
        ELSE   ! scalars
 
+          ib = ifl(lb) 
+          ie = ifu(le)
+          jb = jfl(mb)
+          je = jfu(me)
           DO  l = lb, le
              DO  m = mb, me
@@ -4687,4 +4663,29 @@
 
        ENDIF  ! var
+!
+!--    If the subdomain i- and/or j-dimension (nx/npex and/or ny/npey) is 
+!--    not integer divisible by the grid spacing ratio in its direction,
+!--    the above loops will return with unfilled gaps in the initial fields.
+!--    These gaps, if present, are filled here.  
+       IF  ( ib > ifirst )  THEN
+          DO  i = ifirst, ib-1
+             f(:,:,i) = f(:,:,ib)
+          ENDDO
+       ENDIF
+       IF  ( ie < ilast )  THEN
+          DO  i = ie+1, ilast
+             f(:,:,i) = f(:,:,ie)
+          ENDDO
+       ENDIF
+       IF  ( jb > jfirst )  THEN
+          DO  j = jfirst, jb-1
+             f(:,j,:) = f(:,jb,:)
+          ENDDO
+       ENDIF
+       IF  ( je < jlast )  THEN
+          DO  j = je+1, jlast
+             f(:,j,:) = f(:,je,:)
+          ENDDO
+       ENDIF
 
     END SUBROUTINE pmci_interp_1sto_all
@@ -5754,7 +5755,5 @@
 !
 !--   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.
-
+!--   conditions. This subroutine handles the left and right boundaries. 
       IMPLICIT NONE
 
@@ -5865,7 +5864,7 @@
       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
+!      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
 
 ! 
@@ -5940,101 +5939,111 @@
          var_flag = 0
       ENDIF
-
-      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
+!
+!--   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
+         IF  ( bc_dirichlet_s )  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  ( bc_dirichlet_n )  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:jcnw,0:2)                                &
-                 = fc(0:cg%nz+1,jcsw:jcnw,lbeg:lbeg+2)            
+            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
-
-         IF  ( var == 'v' )  THEN
+!
+!--      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) 
-                  
+         DO  m = jcsw, jcnw-1
+            DO n = 0, kct 
+!
+!--            Use averages of the neighbouring matching grid-line values
+               DO  j = jfl(m), jfl(m+1)
+                  f(kfl(n):kfu(n),j,ibc) = 0.5_wp * ( workarrc_lr(n,m,lw)       &
+                       +  workarrc_lr(n,m+1,lw) )
+               ENDDO
+!
+!--            Then set the values along the matching grid-lines  
+               IF  ( MOD( jfl(m), jgsr ) == 0 )  THEN
+                  f(kfl(n):kfu(n),jfl(m),ibc) = workarrc_lr(n,m,lw)
+               ENDIF
+            ENDDO
+         ENDDO
+!
+!--      Finally, set the values along the last matching grid-line  
+         IF  ( MOD( jfl(jcnw), jgsr ) == 0 )  THEN
+            DO  n = 0, kct
+               f(kfl(n):kfu(n),jfl(jcnw),ibc) = workarrc_lr(n,jcnw,lw)
+            ENDDO
+         ENDIF
+!
+!--      A gap may still remain in some cases if the subdomain size is not 
+!--      divisible by the grid-spacing ratio. In such a case, fill the 
+!--      gap. Note however, this operation may produce some additional 
+!--      momentum conservation error.
+         IF  ( jfl(jcnw) < nyn )  THEN
+            DO  n = 0, kct
+               DO  j = jfl(jcnw)+1, nyn
+                  f(kfl(n):kfu(n),j,ibc) = f(kfl(n):kfu(n),jfl(jcnw),ibc)
                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
+         ENDIF
+
+      ELSE IF  ( var == 'w' )  THEN
+
+         DO  m = jcsw, jcnw
+            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, jcnw
+            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
-
-         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
-
-         ENDIF
+         ENDDO
 
       ELSE   ! u or scalars
          
-         DO  m = jcsw+1, jcnw-1
+         DO  m = jcsw, jcnw
             DO n = 0, kct 
                
                DO  j = jfl(m), jfu(m)
                   DO  k = kfl(n), kfu(n)
-                     f(k,j,ibc) = fc(n,m,l)
+                     f(k,j,ibc) = workarrc_lr(n,m,lw)
                   ENDDO
                ENDDO
@@ -6069,6 +6078,4 @@
 !--   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
 
@@ -6172,7 +6179,7 @@
       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
+!      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
 
 !
@@ -6247,101 +6254,112 @@
          var_flag = 0
       ENDIF
-
-      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
+!
+!--   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      
+         IF  ( bc_dirichlet_l )  THEN
+            workarrc_sn(0:cg%nz+1,0:2,iclw:icrw-1)                              &
+                 = fc(0:cg%nz+1,mbeg:mbeg+2,iclw:icrw-1)
+         ELSE IF  ( bc_dirichlet_r )  THEN
+            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
-
-         IF  ( var == 'u' )  THEN
+!
+!--      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                   
+         DO  l = iclw, icrw-1
+            DO n = 0, kct 
+!
+!--            Use averages of the neighbouring matching grid-line values
+               DO  i = ifl(l), ifl(l+1)
+                  f(kfl(n):kfu(n),jbc,i) = 0.5_wp * ( workarrc_sn(n,mw,l)       &
+                       +  workarrc_sn(n,mw,l+1) )
+               ENDDO
+!
+!--            Then set the values along the matching grid-lines  
+               IF  ( MOD( ifl(l), igsr ) == 0 )  THEN
                   f(kfl(n):kfu(n),jbc,ifl(l)) = workarrc_sn(n,mw,l)
-
+               ENDIF
+
+            ENDDO
+         ENDDO
+!
+!--      Finally, set the values along the last matching grid-line  
+         IF  ( MOD( ifl(icrw), igsr ) == 0 )  THEN
+            DO  n = 0, kct
+               f(kfl(n):kfu(n),jbc,ifl(icrw)) = workarrc_sn(n,mw,icrw)
+            ENDDO
+         ENDIF
+!
+!--      A gap may still remain in some cases if the subdomain size is not 
+!--      divisible by the grid-spacing ratio. In such a case, fill the 
+!--      gap. Note however, this operation may produce some additional 
+!--      momentum conservation error.
+         IF  ( ifl(icrw) < nxr )  THEN
+            DO  n = 0, kct
+               DO  i = ifl(icrw)+1, nxr
+                  f(kfl(n):kfu(n),jbc,i) = f(kfl(n):kfu(n),jbc,ifl(icrw))
                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
-
+         ENDIF
+
+      ELSE IF  ( var == 'w' )  THEN
+
+         DO  l = iclw, icrw
+            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, icrw
+            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
-
-         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
-
-         ENDIF
+         ENDDO
 
       ELSE   ! v or scalars
          
-         DO  l = iclw+1, icrw-1
+         DO  l = iclw, icrw
             DO n = 0, kct 
                
                DO  i = ifl(l), ifu(l)
                   DO  k = kfl(n), kfu(n)
-                     f(k,jbc,i) = fc(n,m,l)
+                     f(k,jbc,i) = workarrc_sn(n,mw,l)
                   ENDDO
                ENDDO
@@ -6374,6 +6392,4 @@
 !--   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
 
@@ -6459,7 +6475,7 @@
       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
+!      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
 
 
@@ -6481,5 +6497,5 @@
          var_flag = 0
       ENDIF
-      n  = kc(k) + noff
+      n  = kc(k) + noff     ! Chance this to get rid of kc. 
       nw = noff 
 !      write(9,"(a1,2x,5(i3,2x))") var, k, kc(k), noff, n, nw 
@@ -6551,10 +6567,92 @@
 #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
-
+      IF  ( var == 'u' )  THEN
+
+         DO  l = iclw, icrw-1
+            DO  m = jcsw, jcnw
+!
+!--            Use averages of the neighbouring matching grid-line values
+               DO  i = ifl(l), ifl(l+1)
+                  f(k,jfl(m):jfu(m),i) = 0.5_wp * ( workarrc_t(nw,m,l)   &
+                       + workarrc_t(nw,m,l+1) )
+               ENDDO
+!
+!--            Then set the values along the matching grid-lines  
+               IF  ( MOD( ifl(l), igsr ) == 0 )  THEN
+                  f(k,jfl(m):jfu(m),ifl(l)) = workarrc_t(nw,m,l)
+               ENDIF
+
+            ENDDO
+         ENDDO
+!
+!--      Finally, set the values along the last matching grid-line  
+         IF  ( MOD( ifl(icrw), igsr ) == 0 )  THEN
+            DO  m = jcsw, jcnw
+               f(k,jfl(m):jfu(m),ifl(icrw)) = workarrc_t(nw,m,icrw)
+            ENDDO
+         ENDIF
+!
+!--      A gap may still remain in some cases if the subdomain size is not 
+!--      divisible by the grid-spacing ratio. In such a case, fill the 
+!--      gap. Note however, this operation may produce some additional 
+!--      momentum conservation error.
+         IF  ( ifl(icrw) < nxr )  THEN
+            DO  m = jcsw, jcnw
+               DO  i = ifl(icrw)+1, nxr
+                  f(k,jfl(m):jfu(m),i) = f(k,jfl(m):jfu(m),ifl(icrw))
+               ENDDO
+            ENDDO
+         ENDIF
+
+      ELSE IF  ( var == 'v' )  THEN
+
+         DO  l = iclw, icrw
+            DO  m = jcsw, jcnw-1
+!
+!--            Use averages of the neighbouring matching grid-line values
+               DO  j = jfl(m), jfl(m+1)                  
+                  f(k,j,ifl(l):ifu(l)) = 0.5_wp * ( workarrc_t(nw,m,l)   &
+                       + workarrc_t(nw,m+1,l) )
+               ENDDO
+!
+!--            Then set the values along the matching grid-lines  
+               IF  ( MOD( jfl(m), jgsr ) == 0 )  THEN
+                  f(k,jfl(m),ifl(l):ifu(l)) = workarrc_t(nw,m,l)
+               ENDIF
+               
+            ENDDO
+
+         ENDDO
+!
+!--      Finally, set the values along the last matching grid-line  
+         IF  ( MOD( jfl(jcnw), jgsr ) == 0 )  THEN
+            DO  l = iclw, icrw
+               f(k,jfl(jcnw),ifl(l):ifu(l)) = workarrc_t(nw,jcnw,l)
+            ENDDO
+         ENDIF
+!
+!--      A gap may still remain in some cases if the subdomain size is not 
+!--      divisible by the grid-spacing ratio. In such a case, fill the 
+!--      gap. Note however, this operation may produce some additional 
+!--      momentum conservation error.
+         IF  ( jfl(jcnw) < nyn )  THEN
+            DO  l = iclw, icrw
+               DO  j = jfl(jcnw)+1, nyn
+                  f(k,j,ifl(l):ifu(l)) = f(k,jfl(jcnw),ifl(l):ifu(l))
+               ENDDO
+            ENDDO
+         ENDIF
+
+      ELSE  ! w or any scalar
+
+         DO  l = iclw, icrw
+            DO  m = jcsw, jcnw
+!
+!--            3rd-order upwind biased interpolation. 
+!               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)
@@ -6565,115 +6663,6 @@
             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) = f_interp
-                  ENDDO
-               ENDDO
-
-            ENDDO
-         ENDDO
-      ENDIF
+
+      ENDIF  ! var
 !
 !--   Just fill up the redundant second ghost-node layer in case of var == w.
