Index: palm/trunk/SOURCE/pmc_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3876)
+++ palm/trunk/SOURCE/pmc_interface_mod.f90	(revision 3883)
@@ -25,4 +25,10 @@
 ! -----------------
 ! $Id$
+! Checks and error messages improved and extended. All the child index bounds in the
+! parent-grid index space are made module variables. Function get_number_of_childs
+! renamed get_number_of_children. A number of variables renamed 
+! and qite a lot of other code reshaping made all around the module. 
+! 
+! 3876 2019-04-08 18:41:49Z knoop
 ! Implemented nesting for salsa variables.
 ! 
@@ -485,18 +491,20 @@
     INTEGER(iwp), PARAMETER ::  child_to_parent = 2   !<
     INTEGER(iwp), PARAMETER ::  parent_to_child = 1   !<
-    INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries (maybe to be made user parameter)
-    INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary (maybe to be made user parameter)
+    INTEGER(iwp), PARAMETER ::  interpolation_scheme_lrsn  = 2  !< Interpolation scheme to be used on lateral boundaries
+    INTEGER(iwp), PARAMETER ::  interpolation_scheme_t     = 3  !< Interpolation scheme to be used on top boundary
 !
 !-- Coupler setup
     INTEGER(iwp), SAVE      ::  comm_world_nesting    !<
     INTEGER(iwp), SAVE      ::  cpl_id  = 1           !<
-    CHARACTER(LEN=32), SAVE ::  cpl_name              !<
     INTEGER(iwp), SAVE      ::  cpl_npe_total         !<
     INTEGER(iwp), SAVE      ::  cpl_parent_id         !<
+    
+    CHARACTER(LEN=32), SAVE ::  cpl_name              !<
+
 !
 !-- Control parameters
+    INTEGER(iwp),     SAVE ::  anterpolation_buffer_width = 2       !< Boundary buffer width for anterpolation
     CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !< steering parameter for data-transfer mode
     CHARACTER(LEN=8), SAVE ::  nesting_mode = 'two-way'             !< steering parameter for 1- or 2-way nesting
-    INTEGER(iwp),     SAVE ::  anterpolation_buffer_width = 2       !< Boundary buffer width for anterpolation
     
     LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
@@ -511,6 +519,4 @@
 !-- 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_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 !< Parent-grid array on child domain - dissipation rate
@@ -529,24 +535,36 @@
     INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
 
-    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c !< coarse grid array on child domain - chemical species
-
-    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_mass_c   !< aerosol mass
-    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_number_c !< aerosol number 
-    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gas_c      !<  salsa gases
+    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c  !< Parent-grid array on child domain - chemical species
+
+    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_mass_c    !< Aerosol mass
+    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  aerosol_number_c  !< Aerosol number 
+    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  salsa_gas_c       !< SALSA gases
 !
 !-- Grid-spacing ratios.
-    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 ::  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
 !
 !-- Global parent-grid index bounds
-    INTEGER(iwp), SAVE ::  iplg     !< Leftmost parent-grid array ip index of the whole child domain
-    INTEGER(iwp), SAVE ::  iprg     !< Rightmost parent-grid array ip index of the whole child domain
-    INTEGER(iwp), SAVE ::  jpsg     !< Southmost parent-grid array jp index of the whole child domain
-    INTEGER(iwp), SAVE ::  jpng     !< Northmost parent-grid array jp index of the whole child domain
-!
-!-- Local parent-grid index bounds (to be moved here from pmci_setup_child)
-!-- EXPLAIN WHY SEVERAL SETS OF PARENT-GRID INDEX BOUNDS ARE NEEDED.
-    
+    INTEGER(iwp), SAVE ::  iplg    !< Leftmost parent-grid array ip index of the whole child domain
+    INTEGER(iwp), SAVE ::  iprg    !< Rightmost parent-grid array ip index of the whole child domain
+    INTEGER(iwp), SAVE ::  jpsg    !< Southmost parent-grid array jp index of the whole child domain
+    INTEGER(iwp), SAVE ::  jpng    !< Northmost parent-grid array jp index of the whole child domain
+!
+!-- Local parent-grid index bounds. Different sets of index bounds are needed for parent-grid arrays (uc, etc),
+!-- for index mapping arrays (iflu, etc) and for work arrays (workarr_lr, etc). This is because these arrays
+!-- have different dimensions depending on the location of the subdomain relative to boundaries and corners.
+    INTEGER(iwp), SAVE ::  ipl     !< Left index limit for children's parent-grid arrays
+    INTEGER(iwp), SAVE ::  ipla    !< Left index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp), SAVE ::  iplw    !< Left index limit for children's parent-grid work arrays
+    INTEGER(iwp), SAVE ::  ipr     !< Right index limit for children's parent-grid arrays
+    INTEGER(iwp), SAVE ::  ipra    !< Right index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp), SAVE ::  iprw    !< Right index limit for children's parent-grid work arrays
+    INTEGER(iwp), SAVE ::  jpn     !< North index limit for children's parent-grid arrays
+    INTEGER(iwp), SAVE ::  jpna    !< North index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp), SAVE ::  jpnw    !< North index limit for children's parent-grid work arrays
+    INTEGER(iwp), SAVE ::  jps     !< South index limit for children's parent-grid arrays
+    INTEGER(iwp), SAVE ::  jpsa    !< South index limit for allocation of index-mapping and other auxiliary arrays
+    INTEGER(iwp), SAVE ::  jpsw    !< South index limit for children's parent-grid work arrays
 !
 !-- Highest prognostic parent-grid k-indices.
@@ -568,21 +586,26 @@
     INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !< child index indicating upper bound of parent grid box on scalar-grid
 !
-!-- Number of fine-grid nodes inside coarse-grid ij-faces
-!-- to be precomputed for anterpolation.
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid boxes contribution to a parent grid box, u-grid
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid boxes contribution to a parent grid box, v-grid
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid boxes contribution to a parent grid box, w-grid
-    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid boxes contribution to a parent grid box, scalar-grid
+!-- Number of child-grid nodes within anterpolation cells to be precomputed for anterpolation.
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_u  !< number of child grid points contributing to a parent grid
+                                                                   !< node in anterpolation, u-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_v  !< number of child grid points contributing to a parent grid
+                                                                   !< node in anterpolation, v-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_w  !< number of child grid points contributing to a parent grid
+                                                                   !< node in anterpolation, w-grid
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  ijkfc_s  !< number of child grid points contributing to a parent grid
+                                                                   !< node in anterpolation, scalar-grid
 !    
 !-- Work arrays for interpolation and user-defined type definitions for horizontal work-array exchange    
-    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
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_lr
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_sn
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: workarr_t
+
+    INTEGER(iwp) :: workarr_lr_exchange_type
+    INTEGER(iwp) :: workarr_sn_exchange_type
+    INTEGER(iwp) :: workarr_t_exchange_type_x
+    INTEGER(iwp) :: workarr_t_exchange_type_y
  
     INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !<
+
     REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !<
     REAL(wp), DIMENSION(2)              ::  zmax_coarse             !<
@@ -607,5 +630,5 @@
     END TYPE parentgrid_def
 
-    TYPE(parentgrid_def), SAVE, PUBLIC     ::  cg                 !< change to pg
+    TYPE(parentgrid_def), SAVE, PUBLIC     ::  pg                 !< Parent-grid information package of type parentgrid_def
 !
 !-- Variables for particle coupling
@@ -617,5 +640,5 @@
        REAL(wp)                            ::  dy                   !<
        REAL(wp)                            ::  dz                   !<
-       REAL(wp)                            ::  lx_coord, lx_coord_b !<
+       REAL(wp)                            ::  lx_coord, lx_coord_b !<   ! split onto separate lines
        REAL(wp)                            ::  rx_coord, rx_coord_b !<
        REAL(wp)                            ::  sy_coord, sy_coord_b !<
@@ -624,8 +647,8 @@
     END TYPE childgrid_def
 
-    TYPE(childgrid_def), SAVE, ALLOCATABLE, DIMENSION(:), PUBLIC :: childgrid !<
-
-    INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: nr_part  !<
-    INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: part_adr !<
+    TYPE(childgrid_def), SAVE, ALLOCATABLE, DIMENSION(:), PUBLIC ::  childgrid  !<
+
+    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  nr_part  !<
+    INTEGER(idp), ALLOCATABLE,DIMENSION(:,:), PUBLIC,TARGET ::  part_adr !<
 
    
@@ -662,7 +685,7 @@
     END INTERFACE
 
-    INTERFACE get_number_of_childs
-       MODULE PROCEDURE get_number_of_childs
-    END  INTERFACE get_number_of_childs
+    INTERFACE get_number_of_children
+       MODULE PROCEDURE get_number_of_children
+    END  INTERFACE get_number_of_children
 
     INTERFACE get_childid
@@ -693,5 +716,5 @@
     PUBLIC pmci_synchronize
     PUBLIC pmci_set_swaplevel
-    PUBLIC get_number_of_childs, get_childid, get_child_edges, get_child_gridspacing
+    PUBLIC get_number_of_children, get_childid, get_child_edges, get_child_gridspacing
 
 
@@ -792,4 +815,5 @@
     INTEGER(iwp) ::  ncpl   !<  number of nest domains
 
+    
 #if defined( __parallel )
     CALL location_message( 'setup the nested model configuration', .FALSE. )
@@ -797,18 +821,19 @@
 !
 !-- Compute absolute coordinates for all models
-    CALL pmci_setup_coordinates
+    CALL pmci_setup_coordinates         ! CONTAIN THIS 
 !
 !-- Determine the number of coupled arrays
-    CALL pmci_num_arrays
+    CALL pmci_num_arrays                ! CONTAIN THIS
 !
 !-- Initialize the child (must be called before pmc_setup_parent)
-    CALL pmci_setup_child
+!  EXTEND THIS COMMENT EXPLAINEIN WHY IT MUST BE CALLED BEFORE    
+    CALL pmci_setup_child               ! CONTAIN THIS
 !
 !-- Initialize PMC parent
-    CALL pmci_setup_parent
+    CALL pmci_setup_parent              ! CONTAIN THIS
 !
 !-- Check for mismatches between settings of master and child variables
 !-- (e.g., all children have to follow the end_time settings of the root master)
-    CALL pmci_check_setting_mismatches
+    CALL pmci_check_setting_mismatches  ! CONTAIN THIS
 !
 !-- Set flag file for combine_plot_fields for precessing the nest output data
@@ -831,51 +856,77 @@
     IMPLICIT NONE
 
-    CHARACTER(LEN=32) ::  myname
-    INTEGER(iwp) ::  child_id         !<
-    INTEGER(iwp) ::  ib = 1           !< running index for aerosol size bins
-    INTEGER(iwp) ::  ic = 1           !< running index for aerosol mass bins
-    INTEGER(iwp) ::  ig = 1           !< running index for salsa gases
-    INTEGER(iwp) ::  ierr             !<
-    INTEGER(iwp) ::  k                !<
-    INTEGER(iwp) ::  m                !<
-    INTEGER(iwp) ::  mid              !<
-    INTEGER(iwp) ::  mm               !<
-    INTEGER(iwp) ::  n = 1            !< running index for chemical species
-    INTEGER(iwp) ::  nest_overlap     !<
-    INTEGER(iwp) ::  nomatch          !<
-    INTEGER(iwp) ::  nx_cl            !<
-    INTEGER(iwp) ::  ny_cl            !<
-    INTEGER(iwp) ::  nz_cl            !<
-
-    INTEGER(iwp), DIMENSION(5) ::  val    !<
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_xl   !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_xr   !<    
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_ys   !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_yn   !<
-    REAL(wp) ::  cl_height        !< 
-    REAL(wp) ::  dx_cl            !<
-    REAL(wp) ::  dy_cl            !<
-    REAL(wp) ::  dz_cl            !<
-    REAL(wp) ::  left_limit       !<
-    REAL(wp) ::  north_limit      !<
-    REAL(wp) ::  right_limit      !<
-    REAL(wp) ::  south_limit      !<
-    REAL(wp) ::  xez              !<
-    REAL(wp) ::  yez              !<
-    REAL(wp), DIMENSION(5) ::  fval                      !< Array for receiving the child-grid spacings etc 
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !<
-
-!
-!   Initialize the pmc parent
+    INTEGER(iwp) ::  child_id           !< Child id-number for the child m
+    INTEGER(iwp) ::  ierr               !< MPI-error code 
+    INTEGER(iwp) ::  kp                 !< Parent-grid index n the z-direction
+    INTEGER(iwp) ::  lb = 1             !< Running index for aerosol size bins
+    INTEGER(iwp) ::  lc = 1             !< Running index for aerosol mass bins
+    INTEGER(iwp) ::  lg = 1             !< Running index for SALSA gases
+    INTEGER(iwp) ::  m                  !< Loop index over all children of the current parent
+    INTEGER(iwp) ::  msib               !< Loop index over all other children than m in case of siblings (parallel children)
+    INTEGER(iwp) ::  n = 1              !< Running index for chemical species
+    INTEGER(iwp) ::  nest_overlap = 0   !< Tag for parallel child-domains' overlap situation (>0 if overlap found)
+    INTEGER(iwp) ::  nomatch = 0        !< Tag for child-domain mismatch situation (>0 if mismatch found)
+    INTEGER(iwp) ::  nx_child           !< Number of child-grid points in the x-direction
+    INTEGER(iwp) ::  ny_child           !< Number of child-grid points in the y-direction
+    INTEGER(iwp) ::  nz_child           !< Number of child-grid points in the z-direction
+    
+    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for receiving the child-grid dimensions from the children
+    
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_left   !< Minimum x-coordinate of the child domain including the ghost
+                                                           !< point layers 
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_x_right  !< Maximum x-coordinate of the child domain including the ghost
+                                                           !< point layers   
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_south  !< Minimum y-coordinate of the child domain including the ghost
+                                                           !< point layers
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_y_north  !< Maximum y-coordinate of the child domain including the ghost
+                                                           !< point layers
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_x  !< Child domain x-coordinate array
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  child_coord_y  !< Child domain y-coordinate array
+    
+    REAL(wp), DIMENSION(5) ::  child_grid_info  !< Array for receiving the child-grid spacings etc from the children
+    
+    REAL(wp) ::  child_height         !< Height of the child domain defined on the child side as zw(nzt+1)
+    REAL(wp) ::  dx_child             !< Child-grid spacing in the x-direction
+    REAL(wp) ::  dy_child             !< Child-grid spacing in the y-direction
+    REAL(wp) ::  dz_child             !< Child-grid spacing in the z-direction
+    REAL(wp) ::  left_limit           !< Left limit for the absolute x-coordinate of the child left boundary 
+    REAL(wp) ::  north_limit          !< North limit for the absolute y-coordinate of the child north boundary
+    REAL(wp) ::  right_limit          !< Right limit for the absolute x-coordinate of the child right boundary
+    REAL(wp) ::  south_limit          !< South limit for the absolute y-coordinate of the child south boundary 
+    REAL(wp) ::  upper_right_coord_x  !< Absolute x-coordinate of the upper right corner of the child domain 
+    REAL(wp) ::  upper_right_coord_y  !< Absolute y-coordinate of the upper right corner of the child domain  
+    REAL(wp) ::  xez                  !< Minimum separation in the x-direction required between the child and
+                                      !< parent boundaries (left or right)
+    REAL(wp) ::  yez                  !< Minimum separation in the y-direction required between the child and
+                                      !< parent boundaries (south or north)
+    CHARACTER(LEN=32) ::  myname      !< String for variable name such as 'u'
+
+    LOGICAL :: m_left_in_msib         !< Logical auxiliary parameter for the overlap test: true if the left border
+                                      !< of the child m is within the x-range of the child msib
+    LOGICAL :: m_right_in_msib        !< Logical auxiliary parameter for the overlap test: true if the right border
+                                      !< of the child m is within the x-range of the child msib 
+    LOGICAL :: msib_left_in_m         !< Logical auxiliary parameter for the overlap test: true if the left border
+                                      !< of the child msib is within the x-range of the child m
+    LOGICAL :: msib_right_in_m        !< Logical auxiliary parameter for the overlap test: true if the right border
+                                      !< of the child msib is within the x-range of the child m
+    LOGICAL :: m_south_in_msib        !< Logical auxiliary parameter for the overlap test: true if the south border
+                                      !< of the child m is within the y-range of the child msib
+    LOGICAL :: m_north_in_msib        !< Logical auxiliary parameter for the overlap test: true if the north border
+                                      !< of the child m is within the y-range of the child msib 
+    LOGICAL :: msib_south_in_m        !< Logical auxiliary parameter for the overlap test: true if the south border
+                                      !< of the child msib is within the y-range of the child m
+    LOGICAL :: msib_north_in_m        !< Logical auxiliary parameter for the overlap test: true if the north border
+                                      !< of the child msib is within the y-range of the child m
+!
+!-- Initialize the current pmc parent
     CALL pmc_parentinit
 !
-!-- Corners of all children of the present parent
-    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN 
-       ALLOCATE( ch_xl(1:SIZE( pmc_parent_for_child ) - 1) )
-       ALLOCATE( ch_xr(1:SIZE( pmc_parent_for_child ) - 1) )
-       ALLOCATE( ch_ys(1:SIZE( pmc_parent_for_child ) - 1) )
-       ALLOCATE( ch_yn(1:SIZE( pmc_parent_for_child ) - 1) )
+!-- Corners of all children of the present parent. Note that 
+!-- SIZE( pmc_parent_for_child ) = 1 if we have no children.
+    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 )  .AND.  myid == 0 )  THEN 
+       ALLOCATE( child_x_left(1:SIZE( pmc_parent_for_child ) - 1) )
+       ALLOCATE( child_x_right(1:SIZE( pmc_parent_for_child ) - 1) )
+       ALLOCATE( child_y_south(1:SIZE( pmc_parent_for_child ) - 1) )
+       ALLOCATE( child_y_north(1:SIZE( pmc_parent_for_child ) - 1) )
     ENDIF
     IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) )  THEN
@@ -883,5 +934,7 @@
     ENDIF
 !
-!-- Get coordinates from all children
+!-- Get coordinates from all children and check that the children match the parent
+!-- domain and each others. Note that SIZE( pmc_parent_for_child ) = 1 
+!-- if we have no children, thence the loop is not executed at all. 
     DO  m = 1, SIZE( pmc_parent_for_child ) - 1
 
@@ -890,35 +943,30 @@
        IF ( myid == 0 )  THEN
 
-          CALL pmc_recv_from_child( child_id, val,  SIZE(val),  0, 123, ierr )
-          CALL pmc_recv_from_child( child_id, fval, SIZE(fval), 0, 124, ierr )
+          CALL pmc_recv_from_child( child_id, child_grid_dim,  SIZE(child_grid_dim), 0, 123, ierr )
+          CALL pmc_recv_from_child( child_id, child_grid_info, SIZE(child_grid_info), 0, 124, ierr )
          
-          nx_cl     = val(1)
-          ny_cl     = val(2)
-          dx_cl     = fval(3)
-          dy_cl     = fval(4)
-          dz_cl     = fval(5)
-          cl_height = fval(1)
-
-          nz_cl = nz
-!
-!--       Find the highest nest level in the parent grid for the reduced z
+          nx_child     = child_grid_dim(1)
+          ny_child     = child_grid_dim(2)
+          dx_child     = child_grid_info(3)
+          dy_child     = child_grid_info(4)
+          dz_child     = child_grid_info(5)
+          child_height = child_grid_info(1)
+!
+!--       Find the highest child-domain level in the parent grid for the reduced z
 !--       transfer
-          DO  k = 1, nz                 
-             IF ( zw(k) > fval(1) )  THEN
-                nz_cl = k
+          DO  kp = 1, nz                 
+             IF ( zw(kp) > child_height )  THEN
+                nz_child = kp
                 EXIT
              ENDIF
           ENDDO
-          zmax_coarse = fval(1:2)
-          cl_height   = fval(1)
+          zmax_coarse  = child_grid_info(1:2)
 !    
 !--       Get absolute coordinates from the child
-          ALLOCATE( cl_coord_x(-nbgp:nx_cl+nbgp) )
-          ALLOCATE( cl_coord_y(-nbgp:ny_cl+nbgp) )
+          ALLOCATE( child_coord_x(-nbgp:nx_child+nbgp) )
+          ALLOCATE( child_coord_y(-nbgp:ny_child+nbgp) )
          
-          CALL pmc_recv_from_child( child_id, cl_coord_x, SIZE( cl_coord_x ),  &
-               0, 11, ierr )
-          CALL pmc_recv_from_child( child_id, cl_coord_y, SIZE( cl_coord_y ),  &
-               0, 12, ierr )
+          CALL pmc_recv_from_child( child_id, child_coord_x, SIZE( child_coord_x ), 0, 11, ierr )
+          CALL pmc_recv_from_child( child_id, child_coord_y, SIZE( child_coord_y ), 0, 12, ierr )
          
           parent_grid_info_real(1) = lower_left_coord_x
@@ -926,34 +974,39 @@
           parent_grid_info_real(3) = dx
           parent_grid_info_real(4) = dy
-          parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx
-          parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy
+
+          upper_right_coord_x      = lower_left_coord_x + ( nx + 1 ) * dx
+          upper_right_coord_y      = lower_left_coord_y + ( ny + 1 ) * dy
+          parent_grid_info_real(5) = upper_right_coord_x
+          parent_grid_info_real(6) = upper_right_coord_y
           parent_grid_info_real(7) = dz(1)
 
           parent_grid_info_int(1)  = nx
           parent_grid_info_int(2)  = ny
-          parent_grid_info_int(3)  = nz_cl
-!
-!--       Check that the child domain matches parent domain. 
-          nomatch = 0
+          parent_grid_info_int(3)  = nz_child
+!
+!--       Check that the child domain matches its parent domain. 
           IF ( nesting_mode == 'vertical' )  THEN
-             right_limit = parent_grid_info_real(5)
-             north_limit = parent_grid_info_real(6)
-             IF ( ( cl_coord_x(nx_cl+1) /= right_limit ) .OR.                  &
-                  ( cl_coord_y(ny_cl+1) /= north_limit ) )  THEN
+!
+!--          In case of vertical nesting, the lateral boundaries must match exactly. 
+             right_limit = upper_right_coord_x
+             north_limit = upper_right_coord_y
+             IF ( ( child_coord_x(nx_child+1) /= right_limit ) .OR.                                 &
+                  ( child_coord_y(ny_child+1) /= north_limit ) )  THEN
                 nomatch = 1
              ENDIF
           ELSE       
 !
-!--       Check that the child domain is completely inside the parent domain. 
+!--          In case of 3-D nesting, check that the child domain is completely 
+!--          inside its parent domain. 
              xez = ( nbgp + 1 ) * dx 
              yez = ( nbgp + 1 ) * dy 
              left_limit  = lower_left_coord_x + xez
-             right_limit = parent_grid_info_real(5) - xez
+             right_limit = upper_right_coord_x - xez
              south_limit = lower_left_coord_y + yez
-             north_limit = parent_grid_info_real(6) - yez
-             IF ( ( cl_coord_x(0) < left_limit )        .OR.                   &
-                  ( cl_coord_x(nx_cl+1) > right_limit ) .OR.                   &
-                  ( cl_coord_y(0) < south_limit )       .OR.                   &
-                  ( cl_coord_y(ny_cl+1) > north_limit ) )  THEN
+             north_limit = upper_right_coord_y - yez
+             IF ( ( child_coord_x(0) < left_limit )  .OR.                                           &
+                  ( child_coord_x(nx_child+1) > right_limit )  .OR.                                 &
+                  ( child_coord_y(0) < south_limit )  .OR.                                          &
+                  ( child_coord_y(ny_child+1) > north_limit ) )  THEN
                 nomatch = 1
              ENDIF
@@ -963,37 +1016,55 @@
 !--       that the top ghost layer of the child grid does not exceed 
 !--       the parent domain top boundary.
-          IF ( cl_height > zw(nz) ) THEN
+          IF ( child_height > zw(nz) ) THEN
              nomatch = 1
           ENDIF
 !
-!--       Check that parallel nest domains, if any, do not overlap.
-          nest_overlap = 0
-          IF ( SIZE( pmc_parent_for_child ) - 1 > 0 )  THEN
-             ch_xl(m) = cl_coord_x(-nbgp)
-             ch_xr(m) = cl_coord_x(nx_cl+nbgp)
-             ch_ys(m) = cl_coord_y(-nbgp)
-             ch_yn(m) = cl_coord_y(ny_cl+nbgp)
-
-             IF ( m > 1 )  THEN
-                DO mm = 1, m - 1
-                   mid = pmc_parent_for_child(mm)
-!
-!--                Check only different nest levels
-                   IF (m_couplers(child_id)%parent_id /= m_couplers(mid)%parent_id)  THEN
-                      IF ( ( ch_xl(m) < ch_xr(mm) .OR.                         &
-                             ch_xr(m) > ch_xl(mm) )  .AND.                     &
-                           ( ch_ys(m) < ch_yn(mm) .OR.                         &
-                             ch_yn(m) > ch_ys(mm) ) )  THEN
-                         nest_overlap = 1
-                      ENDIF
-                   ENDIF
-                ENDDO
-             ENDIF
-          ENDIF
+!--       If parallel child domains (siblings) do exist ( m > 1 ), 
+!--       check that they do not overlap.
+          child_x_left(m)  = child_coord_x(-nbgp)
+          child_x_right(m) = child_coord_x(nx_child+nbgp)
+          child_y_south(m) = child_coord_y(-nbgp)
+          child_y_north(m) = child_coord_y(ny_child+nbgp)
+
+          IF ( nesting_mode /= 'vertical' )  THEN
+!
+!--          Note that the msib-loop is executed only if ( m > 1 ).  
+!--          Also note that the tests have to be made both ways (m vs msib and msib vs m)
+!--          in order to detect all the possible overlap situations.
+             DO  msib = 1, m - 1
+!
+!--             Set some logical auxiliary parameters to simplify the IF-condition.  
+                m_left_in_msib  = ( child_x_left(m)  >= child_x_left(msib) )  .AND.                 &
+                                  ( child_x_left(m)  <= child_x_right(msib) )
+                m_right_in_msib = ( child_x_right(m) >= child_x_left(msib) )  .AND.                 &
+                                  ( child_x_right(m) <= child_x_right(msib) )
+                msib_left_in_m  = ( child_x_left(msib)  >= child_x_left(m) )  .AND.                 &
+                                  ( child_x_left(msib)  <= child_x_right(m) )
+                msib_right_in_m = ( child_x_right(msib) >= child_x_left(m) )  .AND.                 &
+                                  ( child_x_right(msib) <= child_x_right(m) )
+                m_south_in_msib = ( child_y_south(m) >= child_y_south(msib) )  .AND.                &
+                                  ( child_y_south(m) <= child_y_north(msib) )
+                m_north_in_msib = ( child_y_north(m) >= child_y_south(msib) )  .AND.                &
+                                  ( child_y_north(m) <= child_y_north(msib) )
+                msib_south_in_m = ( child_y_south(msib) >= child_y_south(m) )  .AND.                &
+                                  ( child_y_south(msib) <= child_y_north(m) )
+                msib_north_in_m = ( child_y_north(msib) >= child_y_south(m) )  .AND.                &
+                                  ( child_y_north(msib) <= child_y_north(m) )
+                
+                IF ( ( m_left_in_msib  .OR.  m_right_in_msib  .OR.                                  &
+                       msib_left_in_m  .OR.  msib_right_in_m )                                      &
+                     .AND.                                                                          &
+                     ( m_south_in_msib  .OR.  m_north_in_msib  .OR.                                 &
+                       msib_south_in_m  .OR.  msib_north_in_m ) )  THEN
+                     nest_overlap = 1
+                ENDIF
+
+             ENDDO
+          ENDIF          
 
           CALL set_child_edge_coords
 
-          DEALLOCATE( cl_coord_x )
-          DEALLOCATE( cl_coord_y )
+          DEALLOCATE( child_coord_x )
+          DEALLOCATE( child_coord_y )
 !
 !--       Send information about operating mode (LES or RANS) to child. This will be 
@@ -1002,42 +1073,40 @@
 !
 !--       Send coarse grid information to child
-          CALL pmc_send_to_child( child_id, parent_grid_info_real,             &
-                                   SIZE( parent_grid_info_real ), 0, 21,       &
-                                   ierr )
-          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,       &
-                                   22, ierr )
+          CALL pmc_send_to_child( child_id, parent_grid_info_real,                                  &
+                                  SIZE( parent_grid_info_real ), 0, 21,                             &
+                                  ierr )
+          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,                            &
+                                  22, ierr )
 !
 !--       Send local grid to child
-          CALL pmc_send_to_child( child_id, coord_x, nx+1+2*nbgp, 0, 24,       &
-                                   ierr )
-          CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,       &
-                                   ierr )
+          CALL pmc_send_to_child( child_id, coord_x, nx+1+2*nbgp, 0, 24,                            &
+                                  ierr )
+          CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,                            &
+                                  ierr )
 !
 !--       Also send the dzu-, dzw-, zu- and zw-arrays here
-          CALL pmc_send_to_child( child_id, dzu, nz_cl+1, 0, 26, ierr )
-          CALL pmc_send_to_child( child_id, dzw, nz_cl+1, 0, 27, ierr )
-          CALL pmc_send_to_child( child_id, zu,  nz_cl+2, 0, 28, ierr )
-          CALL pmc_send_to_child( child_id, zw,  nz_cl+2, 0, 29, ierr )
-
-       ENDIF
-
-       CALL MPI_BCAST( nomatch, 1, MPI_INTEGER, 0, comm2d, ierr )
-       IF ( nomatch /= 0 )  THEN
-          WRITE ( message_string, * )  'nested child domain does ',            &
-                                       'not fit into its parent domain'
-          CALL message( 'pmci_setup_parent', 'PA0425', 3, 2, 0, 6, 0 )
-       ENDIF
+          CALL pmc_send_to_child( child_id, dzu, nz_child + 1, 0, 26, ierr )
+          CALL pmc_send_to_child( child_id, dzw, nz_child + 1, 0, 27, ierr )
+          CALL pmc_send_to_child( child_id, zu,  nz_child + 2, 0, 28, ierr )
+          CALL pmc_send_to_child( child_id, zw,  nz_child + 2, 0, 29, ierr )
+
+          IF ( nomatch /= 0 )  THEN
+             WRITE ( message_string, * ) 'nested child domain does not fit into its parent domain'
+             CALL message( 'pmci_setup_parent', 'PA0425', 3, 2, 0, 6, 0 )
+          ENDIF
   
-       CALL MPI_BCAST( nest_overlap, 1, MPI_INTEGER, 0, comm2d, ierr )
-       IF ( nest_overlap /= 0  .AND.  nesting_mode /= 'vertical' )  THEN
-          WRITE ( message_string, * )  'nested parallel child domains overlap'
-          CALL message( 'pmci_setup_parent', 'PA0426', 3, 2, 0, 6, 0 )
-       ENDIF
-      
-       CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
+          IF ( nest_overlap /= 0  .AND.  nesting_mode /= 'vertical' )  THEN
+             WRITE ( message_string, * ) 'nested parallel child domains overlap'
+             CALL message( 'pmci_setup_parent', 'PA0426', 3, 2, 0, 6, 0 )
+          ENDIF
+          
+       ENDIF  ! ( myid == 0 ) 
+
+       CALL MPI_BCAST( nz_child, 1, MPI_INTEGER, 0, comm2d, ierr )
 
        CALL MPI_BCAST( childgrid(m), STORAGE_SIZE(childgrid(1))/8, MPI_BYTE, 0, comm2d, ierr )
 !
 !--    TO_DO: Klaus: please give a comment what is done here
+!  DO IT YOURSELF       
        CALL pmci_create_index_list
 !
@@ -1047,39 +1116,40 @@
 !--    have to be specified again
        CALL pmc_s_clear_next_array_list
-       DO  WHILE ( pmc_s_getnextarray( child_id, myname ) )
+       DO WHILE ( pmc_s_getnextarray( child_id, myname ) )
           IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
              CALL pmci_set_array_pointer( myname, child_id = child_id,         &
-                                          nz_cl = nz_cl, n = n )
+                                          nz_child = nz_child, n = n )
              n = n + 1  
           ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
              CALL pmci_set_array_pointer( myname, child_id = child_id,         &
-                                          nz_cl = nz_cl, n = ib )
-             ib = ib + 1 
+                                          nz_child = nz_child, n = lb )
+             lb = lb + 1 
           ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
              CALL pmci_set_array_pointer( myname, child_id = child_id,         &
-                                          nz_cl = nz_cl, n = ic )
-             ic = ic + 1 
+                                          nz_child = nz_child, n = lc )
+             lc = lc + 1 
           ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.                  &
              .NOT. salsa_gases_from_chem )                                     &
           THEN
              CALL pmci_set_array_pointer( myname, child_id = child_id,         &
-                                          nz_cl = nz_cl, n = ig )
-             ig = ig + 1
+                                          nz_child = nz_child, n = lg )
+             lg = lg + 1
           ELSE
-             CALL pmci_set_array_pointer( myname, child_id = child_id,         &
-                                          nz_cl = nz_cl )
+             CALL pmci_set_array_pointer( myname, child_id = child_id, nz_child = nz_child )
           ENDIF
        ENDDO
 
        CALL pmc_s_setind_and_allocmem( child_id )
-    ENDDO
-
-    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN 
-       DEALLOCATE( ch_xl )
-       DEALLOCATE( ch_xr )
-       DEALLOCATE( ch_ys )
-       DEALLOCATE( ch_yn )
+       
+    ENDDO  ! m
+
+    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN
+       DEALLOCATE( child_x_left )
+       DEALLOCATE( child_x_right )
+       DEALLOCATE( child_y_south )
+       DEALLOCATE( child_y_north )
     ENDIF
 
+    
  CONTAINS
 
@@ -1108,4 +1178,5 @@
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list         !<
 
+       
        IF ( myid == 0 )  THEN
 !          
@@ -1138,8 +1209,8 @@
 !
 !--          Area along y required by actual child PE
-             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)  !: j = jcs, jcn of PE# k 
+             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)  !: j = jps, jpn of PE# k 
 !
 !--             Area along x required by actual child PE
-                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)  !: i = icl, icr of PE# k 
+                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)  !: i = ipl, ipr of PE# k 
 
                    px = i / nrx
@@ -1192,23 +1263,24 @@
         IMPLICIT  NONE
 
-        INTEGER(iwp) :: nbgp_lpm = 1
-
-        nbgp_lpm = min(nbgp_lpm, nbgp)
-
-        childgrid(m)%nx = nx_cl
-        childgrid(m)%ny = ny_cl
-        childgrid(m)%nz = nz_cl
-        childgrid(m)%dx = dx_cl
-        childgrid(m)%dy = dy_cl
-        childgrid(m)%dz = dz_cl
-
-        childgrid(m)%lx_coord   = cl_coord_x(0)
-        childgrid(m)%lx_coord_b = cl_coord_x(-nbgp_lpm)
-        childgrid(m)%rx_coord   = cl_coord_x(nx_cl)+dx_cl
-        childgrid(m)%rx_coord_b = cl_coord_x(nx_cl+nbgp_lpm)+dx_cl
-        childgrid(m)%sy_coord   = cl_coord_y(0)
-        childgrid(m)%sy_coord_b = cl_coord_y(-nbgp_lpm)
-        childgrid(m)%ny_coord   = cl_coord_y(ny_cl)+dy_cl
-        childgrid(m)%ny_coord_b = cl_coord_y(ny_cl+nbgp_lpm)+dy_cl
+        INTEGER(iwp) ::  nbgp_lpm = 1  !<
+
+        
+        nbgp_lpm = MIN( nbgp_lpm, nbgp )
+
+        childgrid(m)%nx = nx_child
+        childgrid(m)%ny = ny_child
+        childgrid(m)%nz = nz_child
+        childgrid(m)%dx = dx_child
+        childgrid(m)%dy = dy_child
+        childgrid(m)%dz = dz_child
+
+        childgrid(m)%lx_coord   = child_coord_x(0)
+        childgrid(m)%lx_coord_b = child_coord_x(-nbgp_lpm)
+        childgrid(m)%rx_coord   = child_coord_x(nx_child) + dx_child
+        childgrid(m)%rx_coord_b = child_coord_x(nx_child+nbgp_lpm) + dx_child
+        childgrid(m)%sy_coord   = child_coord_y(0)
+        childgrid(m)%sy_coord_b = child_coord_y(-nbgp_lpm)
+        childgrid(m)%ny_coord   = child_coord_y(ny_child) + dy_child
+        childgrid(m)%ny_coord_b = child_coord_y(ny_child+nbgp_lpm) + dy_child
         childgrid(m)%uz_coord   = zmax_coarse(2)
         childgrid(m)%uz_coord_b = zmax_coarse(1)
@@ -1226,37 +1298,22 @@
     IMPLICIT NONE
 
-
-    CHARACTER(LEN=da_namelen) ::  myname     !<
-    CHARACTER(LEN=5) ::  salsa_char          !<
-    INTEGER(iwp) ::  ib         !< running index for aerosol size bins
-    INTEGER(iwp) ::  ic         !< running index for aerosol mass bins
-    INTEGER(iwp) ::  ierr       !<
-    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) ::  ig         !< running index for salsa gases
-    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(3) ::  val  !< Array for sending the child-grid dimensions to parent 
-    REAL(wp) ::  xcs        !<
-    REAL(wp) ::  xce        !<
-    REAL(wp) ::  ycs        !<
-    REAL(wp) ::  yce        !<
-
-    REAL(wp), DIMENSION(5) ::  fval     !< Array for sending the child-grid spacings etc to parent 
-                                             
+    INTEGER(iwp) ::  ierr                          !< MPI error code
+    INTEGER(iwp) ::  lb                            !< Running index for aerosol size bins
+    INTEGER(iwp) ::  lc                            !< Running index for aerosol mass bins
+    INTEGER(iwp) ::  lg                            !< Running index for salsa gases
+    INTEGER(iwp) ::  n                             !< Running index for number of chemical species
+    INTEGER(iwp), DIMENSION(3) ::  child_grid_dim  !< Array for sending the child-grid dimensions to parent 
+
+    REAL(wp), DIMENSION(5) ::  child_grid_info     !< Array for sending the child-grid spacings etc to parent 
+         
+    CHARACTER( LEN=da_namelen ) ::  myname         !< 
+    CHARACTER(LEN=5) ::  salsa_char                !<
+    
 !
 !-- Child setup
 !-- Root model does not have a parent and is not a child, therefore no child setup on root model
     IF ( .NOT. pmc_is_rootmodel() )  THEN
-
+!
+!--    ADD A DESCRIPTION HERE WHAT PMC_CHILDINIT DOES       
        CALL pmc_childinit
 !
@@ -1267,5 +1324,5 @@
 !--    in subroutines:
 !--    pmci_set_array_pointer (for parent arrays)
-!--    pmci_create_child_arrays (for child arrays)
+!--    pmci_create_childs_parent_grid_arrays (for child's parent-grid arrays)
        CALL pmc_set_dataarray_name( 'coarse', 'u'  ,'fine', 'u',  ierr )
        CALL pmc_set_dataarray_name( 'coarse', 'v'  ,'fine', 'v',  ierr )
@@ -1307,5 +1364,5 @@
        ENDIF
 
-       IF( particle_advection )  THEN
+       IF ( particle_advection )  THEN
           CALL pmc_set_dataarray_name( 'coarse', 'nr_part'  ,'fine',           &
                'nr_part',  ierr )
@@ -1315,9 +1372,9 @@
        
        IF ( air_chemistry  .AND.  nest_chemistry )  THEN
-          DO  n = 1, nspec
+          DO n = 1, nspec
              CALL pmc_set_dataarray_name( 'coarse',                            &
                                           'chem_' //                           &
                                           TRIM( chem_species(n)%name ),        &
-                                          'fine',                              &
+                                         'fine',                               &
                                           'chem_' //                           &
                                           TRIM( chem_species(n)%name ),        &
@@ -1327,6 +1384,6 @@
 
        IF ( salsa  .AND.  nest_salsa )  THEN
-          DO  ib = 1, nbins_aerosol
-             WRITE(salsa_char,'(i0)') ib
+          DO  lb = 1, nbins_aerosol
+             WRITE(salsa_char,'(i0)') lb
              CALL pmc_set_dataarray_name( 'coarse',                            &
                                           'an_' //                             &
@@ -1337,6 +1394,6 @@
                                           ierr )
           ENDDO
-          DO  ic = 1, nbins_aerosol * ncomponents_mass
-             WRITE(salsa_char,'(i0)') ic
+          DO  lc = 1, nbins_aerosol * ncomponents_mass
+             WRITE(salsa_char,'(i0)') lc
              CALL pmc_set_dataarray_name( 'coarse',                            &
                                           'am_' //                             &
@@ -1348,6 +1405,6 @@
           ENDDO
           IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                WRITE(salsa_char,'(i0)') ig
+             DO  lg = 1, ngases_salsa
+                WRITE(salsa_char,'(i0)') lg
                 CALL pmc_set_dataarray_name( 'coarse',                         &
                                              'sg_' //                          &
@@ -1364,17 +1421,17 @@
 !
 !--    Send grid to parent
-       val(1)  = nx
-       val(2)  = ny
-       val(3)  = nz
-       fval(1) = zw(nzt+1)
-       fval(2) = zw(nzt)
-       fval(3) = dx
-       fval(4) = dy
-       fval(5) = dz(1)
+       child_grid_dim(1)  = nx
+       child_grid_dim(2)  = ny
+       child_grid_dim(3)  = nz
+       child_grid_info(1) = zw(nzt+1)
+       child_grid_info(2) = zw(nzt)
+       child_grid_info(3) = dx
+       child_grid_info(4) = dy
+       child_grid_info(5) = dz(1)
 
        IF ( myid == 0 )  THEN
 
-          CALL pmc_send_to_parent( val, SIZE( val ), 0, 123, ierr )
-          CALL pmc_send_to_parent( fval, SIZE( fval ), 0, 124, ierr )
+          CALL pmc_send_to_parent( child_grid_dim, SIZE( child_grid_dim ), 0, 123, ierr )
+          CALL pmc_send_to_parent( child_grid_info, SIZE( child_grid_info ), 0, 124, ierr )
           CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
           CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
@@ -1393,37 +1450,40 @@
        CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr )
 
-       cg%dx = parent_grid_info_real(3)
-       cg%dy = parent_grid_info_real(4)
-       cg%dz = parent_grid_info_real(7)
-       cg%nx = parent_grid_info_int(1)
-       cg%ny = parent_grid_info_int(2)
-       cg%nz = parent_grid_info_int(3)
+       pg%dx = parent_grid_info_real(3)
+       pg%dy = parent_grid_info_real(4)
+       pg%dz = parent_grid_info_real(7)
+       pg%nx = parent_grid_info_int(1)
+       pg%ny = parent_grid_info_int(2)
+       pg%nz = parent_grid_info_int(3)
 !
 !--    Get parent coordinates on coarse grid
-       ALLOCATE( cg%coord_x(-nbgp:cg%nx+nbgp) )
-       ALLOCATE( cg%coord_y(-nbgp:cg%ny+nbgp) )
-       ALLOCATE( cg%dzu(1:cg%nz+1) )
-       ALLOCATE( cg%dzw(1:cg%nz+1) )
-       ALLOCATE( cg%zu(0:cg%nz+1) )
-       ALLOCATE( cg%zw(0:cg%nz+1) )
+       ALLOCATE( pg%coord_x(-nbgp:pg%nx+nbgp) )
+       ALLOCATE( pg%coord_y(-nbgp:pg%ny+nbgp) )
+       ALLOCATE( pg%dzu(1:pg%nz+1) )
+       ALLOCATE( pg%dzw(1:pg%nz+1) )
+       ALLOCATE( pg%zu(0:pg%nz+1) )
+       ALLOCATE( pg%zw(0:pg%nz+1) )
 !
 !--    Get coarse grid coordinates and values of the z-direction from the parent
        IF ( myid == 0)  THEN
-          CALL pmc_recv_from_parent( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr )
-          CALL pmc_recv_from_parent( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr )
-          CALL pmc_recv_from_parent( cg%dzu, cg%nz+1, 0, 26, ierr )
-          CALL pmc_recv_from_parent( cg%dzw, cg%nz+1, 0, 27, ierr )
-          CALL pmc_recv_from_parent( cg%zu, cg%nz+2, 0, 28, ierr )
-          CALL pmc_recv_from_parent( cg%zw, cg%nz+2, 0, 29, ierr )
+          CALL pmc_recv_from_parent( pg%coord_x, pg%nx+1+2*nbgp, 0, 24, ierr )
+          CALL pmc_recv_from_parent( pg%coord_y, pg%ny+1+2*nbgp, 0, 25, ierr )
+          CALL pmc_recv_from_parent( pg%dzu, pg%nz+1, 0, 26, ierr )
+          CALL pmc_recv_from_parent( pg%dzw, pg%nz+1, 0, 27, ierr )
+          CALL pmc_recv_from_parent( pg%zu, pg%nz+2, 0, 28, ierr )
+          CALL pmc_recv_from_parent( pg%zw, pg%nz+2, 0, 29, ierr )
        ENDIF
 !
 !--    Broadcast this information
-       CALL MPI_BCAST( cg%coord_x, cg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( cg%coord_y, cg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( cg%dzu, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( cg%dzw, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( cg%zu, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
-       CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%coord_x, pg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%coord_y, pg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%dzu, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%dzw, pg%nz+1, MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%zu, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
+       CALL MPI_BCAST( pg%zw, pg%nz+2,  MPI_REAL, 0, comm2d, ierr )
        CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )
+
+!  CHECK IF pmci_check_grid_matching COULD BE MOVED HERE.
+       
 !
 !--    Find the index bounds for the nest domain in the coarse-grid index space
@@ -1437,8 +1497,8 @@
        CALL  pmc_c_clear_next_array_list
 
-       ib = 1
-       ic = 1
-       ig = 1
        n  = 1
+       lb = 1
+       lc = 1
+       lg = 1
 
        DO  WHILE ( pmc_c_getnextarray( myname ) )
@@ -1451,22 +1511,17 @@
 !--       species and increment this subsequently.
           IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
-             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz, n )
+             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, n )
              n = n + 1    
           ELSEIF ( INDEX( TRIM( myname ), 'an_' ) /= 0 )  THEN
-             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz,&
-                                             ib )
-             ib = ib + 1
+             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lb )
+             lb = lb + 1
           ELSEIF ( INDEX( TRIM( myname ), 'am_' ) /= 0 )  THEN
-             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz,&
-                                             ic )
-             ic = ic + 1
-          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.                  &
-             .NOT. salsa_gases_from_chem )                                     &
-          THEN
-             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz,&
-                                             ig )
-             ig = ig + 1
+             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lc )
+             lc = lc + 1
+          ELSEIF ( INDEX( TRIM( myname ), 'sg_' ) /= 0  .AND.  .NOT.  salsa_gases_from_chem )  THEN
+             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz, lg )
+             lg = lg + 1
           ELSE
-             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz )
+             CALL pmci_create_childs_parent_grid_arrays ( myname, ipl, ipr, jps, jpn, pg%nz )
           ENDIF
        ENDDO
@@ -1491,17 +1546,24 @@
 
        INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all     !< Transfer array for parent-grid index bounds
+
        INTEGER(iwp), DIMENSION(4)          ::  parent_bound_global  !< Transfer array for global parent-grid index bounds
        INTEGER(iwp), DIMENSION(2)          ::  size_of_array        !<
-       INTEGER(iwp) :: i        !< 
-       INTEGER(iwp) :: iauxl    !<
-       INTEGER(iwp) :: iauxr    !<
-       INTEGER(iwp) :: ijaux    !< 
-       INTEGER(iwp) :: j        !<
-       INTEGER(iwp) :: jauxs    !<
-       INTEGER(iwp) :: jauxn    !< 
+
+       INTEGER(iwp) ::  i       !< 
+       INTEGER(iwp) ::  iauxl   !<
+       INTEGER(iwp) ::  iauxr   !<
+       INTEGER(iwp) ::  ijaux   !< 
+       INTEGER(iwp) ::  j       !<
+       INTEGER(iwp) ::  jauxs   !<
+       INTEGER(iwp) ::  jauxn   !< 
+
        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
+       REAL(wp) ::  xcs         !< RENAME 
+       REAL(wp) ::  xce         !< RENAME 
+       REAL(wp) ::  ycs         !< RENAME 
+       REAL(wp) ::  yce         !< RENAME 
 
 !
@@ -1511,9 +1573,11 @@
 !--    Else the parent-grid cell is included in the neighbouring subdomain's 
 !--    anterpolation domain, or not included at all if we are at the outer 
-!--    edge of the child domain.
+!--    edge of the child domain. This may occur especially when a large grid-spacing
+!--    ratio is used.       
 !
 !--    Left
-       IF  ( bc_dirichlet_l )  THEN
-          xexl  = 2 * cg%dx
+!--    EXPLAIN THE EXTENSION HERE AND IN THE OTHER OCCASIONS (r, s, n)       
+       IF ( bc_dirichlet_l )  THEN
+          xexl  = 2 * pg%dx
           iauxl = 0
        ELSE
@@ -1522,7 +1586,7 @@
        ENDIF
        xcs     = coord_x(nxl) - xexl
-       DO  i = 0, cg%nx
-          IF ( cg%coord_x(i) + 0.5_wp * cg%dx >= xcs )  THEN
-             icl = MAX( 0, i )
+       DO  i = 0, pg%nx
+          IF ( pg%coord_x(i) + 0.5_wp * pg%dx >= xcs )  THEN   ! Consider changing >= to ==
+             ipl = MAX( 0, i )
              EXIT
           ENDIF
@@ -1530,6 +1594,6 @@
 !
 !--    Right
-       IF  ( bc_dirichlet_r )  THEN
-          xexr  = 2 * cg%dx
+       IF ( bc_dirichlet_r )  THEN
+          xexr  = 2 * pg%dx
           iauxr = 0  
        ELSE
@@ -1538,7 +1602,7 @@
        ENDIF
        xce  = coord_x(nxr+1) + xexr
-       DO  i = cg%nx, 0 , -1
-          IF ( cg%coord_x(i) + 0.5_wp * cg%dx <= xce )  THEN
-             icr = MIN( cg%nx, MAX( icl, i ) )
+       DO  i = pg%nx, 0 , -1
+          IF ( pg%coord_x(i) + 0.5_wp * pg%dx <= xce )  THEN
+             ipr = MIN( pg%nx, MAX( ipl, i ) )
              EXIT
           ENDIF
@@ -1546,6 +1610,6 @@
 !
 !--    South
-       IF  ( bc_dirichlet_s )  THEN
-          yexs  = 2 * cg%dy
+       IF ( bc_dirichlet_s )  THEN
+          yexs  = 2 * pg%dy
           jauxs = 0  
        ELSE
@@ -1554,7 +1618,7 @@
        ENDIF
        ycs  = coord_y(nys) - yexs
-       DO  j = 0, cg%ny
-          IF ( cg%coord_y(j) + 0.5_wp * cg%dy >= ycs )  THEN
-             jcs = MAX( 0, j )
+       DO  j = 0, pg%ny
+          IF ( pg%coord_y(j) + 0.5_wp * pg%dy >= ycs )  THEN
+             jps = MAX( 0, j )
              EXIT
           ENDIF
@@ -1563,5 +1627,5 @@
 !--    North
        IF  ( bc_dirichlet_n )  THEN
-          yexn  = 2 * cg%dy
+          yexn  = 2 * pg%dy
           jauxn = 0
        ELSE
@@ -1570,53 +1634,55 @@
        ENDIF
        yce  = coord_y(nyn+1) + yexn
-       DO  j = cg%ny, 0 , -1
-          IF ( cg%coord_y(j) + 0.5_wp * cg%dy <= yce )  THEN
-             jcn = MIN( cg%ny, MAX( jcs, j ) )
+       DO  j = pg%ny, 0 , -1
+          IF ( pg%coord_y(j) + 0.5_wp * pg%dy <= yce )  THEN
+             jpn = MIN( pg%ny, MAX( jps, j ) )
              EXIT
           ENDIF
        ENDDO
 !
-!--    Make sure that the indexing is contiguous (no gaps, no overlaps) 
-#if defined( __parallel )
+!--    Make sure that the indexing is contiguous (no gaps, no overlaps). 
+!--    This is a safety measure mainly for cases with high grid-spacing 
+!--    ratio and narrow child subdomains.
        IF ( nxl == 0 )  THEN
-          CALL MPI_SEND( icr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+          CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
        ELSE IF ( nxr == nx )  THEN
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr )
-          icl = ijaux + 1
+          ipl = ijaux + 1
        ELSE
-          CALL MPI_SEND( icr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
+          CALL MPI_SEND( ipr, 1, MPI_INTEGER, pright, 717, comm2d, ierr )
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, pleft, 717, comm2d, status, ierr ) 
-          icl = ijaux + 1
+          ipl = ijaux + 1
        ENDIF
        IF ( nys == 0 )  THEN
-          CALL MPI_SEND( jcn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+          CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
        ELSE IF ( nyn == ny )  THEN
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr )
-          jcs = ijaux + 1
+          jps = ijaux + 1
        ELSE
-          CALL MPI_SEND( jcn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
+          CALL MPI_SEND( jpn, 1, MPI_INTEGER, pnorth, 719, comm2d, ierr )
           CALL MPI_RECV( ijaux, 1, MPI_INTEGER, psouth, 719, comm2d, status, ierr ) 
-          jcs = ijaux + 1
-       ENDIF
-#endif       
-
-       WRITE(9,"('Pmci_map_fine_to_coarse_grid. Parent-grid array bounds: ',4(i4,2x))") icl, icr, jcs, jcn
+          jps = ijaux + 1
+       ENDIF
+
+       WRITE(9,"('Pmci_map_fine_to_coarse_grid. Parent-grid array bounds: ',4(i4,2x))")             &
+            ipl, ipr, jps, jpn
        FLUSH(9)
 
-       coarse_bound(1) = icl
-       coarse_bound(2) = icr
-       coarse_bound(3) = jcs
-       coarse_bound(4) = jcn
+       coarse_bound(1) = ipl
+       coarse_bound(2) = ipr
+       coarse_bound(3) = jps
+       coarse_bound(4) = jpn
        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
+       ipla = ipl - iauxl
+       ipra = ipr + iauxr
+       jpsa = jps - jauxs
+       jpna = jpn + jauxn
 !
 !--    Note that MPI_Gather receives data from all processes in the rank order
 !--    This fact is exploited in creating the index list in pmci_create_index_list
-       CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5,     &
+!    IMPROVE THIS COMMENT. EXPLAIN WHERE THIS INFORMATION IS NEEDED.
+       CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5,                          &
                         MPI_INTEGER, 0, comm2d, ierr )
 
@@ -1625,6 +1691,5 @@
           size_of_array(2) = SIZE( coarse_bound_all, 2 )
           CALL pmc_send_to_parent( size_of_array, 2, 0, 40, ierr )
-          CALL pmc_send_to_parent( coarse_bound_all, SIZE( coarse_bound_all ), &
-               0, 41, ierr )
+          CALL pmc_send_to_parent( coarse_bound_all, SIZE( coarse_bound_all ), 0, 41, ierr )
 !
 !--       Determine the global parent-grid index bounds       
@@ -1641,6 +1706,7 @@
        jpsg = parent_bound_global(3)
        jpng = parent_bound_global(4)
-       WRITE(9,"('Pmci_map_fine_to_coarse_grid. Global parent-grid index bounds iplg, iprg, jpsg, jpng: ',4(i4,2x))") iplg, iprg, jpsg, jpng
-       FLUSH(9)
+       WRITE( 9, "('Pmci_map_fine_to_coarse_grid. Global parent-grid index bounds iplg, iprg, jpsg, jpng: ',4(i4,2x))" ) &
+            iplg, iprg, jpsg, jpng
+       FLUSH( 9 )
        
     END SUBROUTINE pmci_map_fine_to_coarse_grid
@@ -1654,24 +1720,24 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i         !< Child-grid index
-       INTEGER(iwp) ::  ii        !< Parent-grid index
+       INTEGER(iwp) ::  i         !< Child-grid index in the x-direction
+       INTEGER(iwp) ::  ii        !< Parent-grid index in the x-direction
        INTEGER(iwp) ::  istart    !<
        INTEGER(iwp) ::  ir        !<
-       INTEGER(iwp) ::  iw        !< Child-grid index limited to -1 <= iw <= nx+1
-       INTEGER(iwp) ::  j         !< Child-grid index
-       INTEGER(iwp) ::  jj        !< Parent-grid index
+       INTEGER(iwp) ::  iw        !< Child-grid index limited to -1 <= iw <= nx+1 for wall_flags_0
+       INTEGER(iwp) ::  j         !< Child-grid index in the y-direction
+       INTEGER(iwp) ::  jj        !< Parent-grid index in the y-direction
        INTEGER(iwp) ::  jstart    !<
        INTEGER(iwp) ::  jr        !<
-       INTEGER(iwp) ::  jw        !< Child-grid index limited to -1 <= jw <= ny+1
-       INTEGER(iwp) ::  k         !< Child-grid index
-       INTEGER(iwp) ::  kk        !< Parent-grid index
+       INTEGER(iwp) ::  jw        !< Child-grid index limited to -1 <= jw <= ny+1 for wall_flags_0
+       INTEGER(iwp) ::  k         !< Child-grid index in the z-direction
+       INTEGER(iwp) ::  kk        !< Parent-grid index in the z-direction
        INTEGER(iwp) ::  kstart    !<
-       INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1
+       INTEGER(iwp) ::  kw        !< Child-grid index limited to kw <= nzt+1 for wall_flags_0
      
 !
 !--    Allocate child-grid work arrays for interpolation.
-       igsr = NINT( cg%dx / dx, iwp )
-       jgsr = NINT( cg%dy / dy, iwp )
-       kgsr = NINT( cg%dzw(1) / dzw(1), iwp )
+       igsr = NINT( pg%dx / dx, iwp )
+       jgsr = NINT( pg%dy / dy, iwp )
+       kgsr = NINT( pg%dzw(1) / dzw(1), iwp )
        WRITE(9,"('igsr, jgsr, kgsr: ',3(i3,2x))") igsr, jgsr, kgsr
        FLUSH(9)
@@ -1688,5 +1754,5 @@
 !--    coarse-grid levels below the child top-boundary level.
        kk = 0
-       DO  WHILE ( cg%zu(kk) <= zu(nzt) )
+       DO WHILE ( pg%zu(kk) <= zu(nzt) )
           kk = kk + 1
        ENDDO
@@ -1694,32 +1760,28 @@
 
        kk = 0
-       DO  WHILE ( cg%zw(kk) <= zw(nzt-1) )
+       DO WHILE ( pg%zw(kk) <= zw(nzt-1) )
           kk = kk + 1
        ENDDO
        kctw = kk - 1
 
-       WRITE(9,"('kcto, kctw = ', 2(i3,2x))") kcto, kctw
-       FLUSH(9)
+       WRITE( 9, "('kcto, kctw = ', 2(i3,2x))" ) kcto, kctw
+       FLUSH( 9 )
        
-       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) )       
-       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,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) )
+       ALLOCATE( iflu(ipla:ipra) )
+       ALLOCATE( iflo(ipla:ipra) )
+       ALLOCATE( ifuu(ipla:ipra) )
+       ALLOCATE( ifuo(ipla:ipra) )
+       ALLOCATE( jflv(jpsa:jpna) )
+       ALLOCATE( jflo(jpsa:jpna) )
+       ALLOCATE( jfuv(jpsa:jpna) )
+       ALLOCATE( jfuo(jpsa:jpna) )       
+       ALLOCATE( kflw(0:pg%nz+1) )
+       ALLOCATE( kflo(0:pg%nz+1) )
+       ALLOCATE( kfuw(0:pg%nz+1) )
+       ALLOCATE( kfuo(0:pg%nz+1) )
+       ALLOCATE( ijkfc_u(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
+       ALLOCATE( ijkfc_v(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
+       ALLOCATE( ijkfc_w(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
+       ALLOCATE( ijkfc_s(0:pg%nz+1,jpsa:jpna,ipla:ipra) )
 
        ijkfc_u = 0
@@ -1730,10 +1792,13 @@
 !--    i-indices of u for each ii-index value
        istart = nxlg
-       DO  ii = icla, icra
+       DO  ii = ipla, ipra
 !
 !--       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.
+!--       However, icru still has to be stored separately as these index bounds
+!--       are passed as arguments to the interpolation and anterpolation 
+!--       subroutines.
           i = istart
-          DO WHILE ( coord_x(i) < cg%coord_x(ii) .AND. i < nxrg )
+          DO WHILE ( coord_x(i) < pg%coord_x(ii) .AND. i < nxrg )
              i = i + 1
           ENDDO
@@ -1743,15 +1808,15 @@
 !
 !--       Print out the index bounds for checking and debugging purposes
-          WRITE(9,"('pmci_define_index_mapping, ii, iflu, ifuu: ', 3(i4,2x))")    &
+          WRITE( 9, "('pmci_define_index_mapping, ii, iflu, ifuu: ', 3(i4,2x))" )                   &
                ii, iflu(ii), ifuu(ii)
-          FLUSH(9)
+          FLUSH( 9 )
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    i-indices of others for each ii-index value
        istart = nxlg
-       DO  ii = icla, icra
+       DO  ii = ipla, ipra
           i = istart
-          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) )  .AND.     &
+          DO WHILE ( ( coord_x(i) + 0.5_wp * dx < pg%coord_x(ii) )  .AND.     &
                       ( i < nxrg ) )
              i  = i + 1
@@ -1759,5 +1824,5 @@
           iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
           ir = i
-          DO  WHILE ( ( coord_x(ir) + 0.5_wp * dx <= cg%coord_x(ii) + cg%dx )  &
+          DO WHILE ( ( coord_x(ir) + 0.5_wp * dx <= pg%coord_x(ii) + pg%dx )  &
                       .AND.  ( i < nxrg+1 ) )
              i  = i + 1
@@ -1768,18 +1833,21 @@
 !
 !--       Print out the index bounds for checking and debugging purposes
-          WRITE(9,"('pmci_define_index_mapping, ii, iflo, ifuo: ', 3(i4,2x))")    &
+          WRITE( 9, "('pmci_define_index_mapping, ii, iflo, ifuo: ', 3(i4,2x))" )                   &
                ii, iflo(ii), ifuo(ii)
-          FLUSH(9)
+          FLUSH( 9 )
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    j-indices of v for each jj-index value
        jstart = nysg
-       DO  jj = jcsa, jcna
+       DO  jj = jpsa, jpna
 !
 !--       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.
+!--       However, jcnv still has to be stored separately as these index bounds 
+!--       are passed as arguments to the interpolation and anterpolation 
+!--       subroutines.
           j = jstart
-          DO WHILE ( coord_y(j) < cg%coord_y(jj) .AND. j < nyng )
+          DO WHILE ( coord_y(j) < pg%coord_y(jj) .AND. j < nyng )
              j = j + 1
           ENDDO
@@ -1789,22 +1857,20 @@
 !
 !--       Print out the index bounds for checking and debugging purposes
-          WRITE(9,"('pmci_define_index_mapping, jj, jflv, jfuv: ', 3(i4,2x))")    &
+          WRITE( 9, "('pmci_define_index_mapping, jj, jflv, jfuv: ', 3(i4,2x))" )                   &
                jj, jflv(jj), jfuv(jj)
           FLUSH(9)
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    j-indices of others for each jj-index value
        jstart = nysg
-       DO  jj = jcsa, jcna
+       DO  jj = jpsa, jpna
           j = jstart
-          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) )  .AND.     &
-                      ( j < nyng ) )
+          DO WHILE ( ( coord_y(j) + 0.5_wp * dy < pg%coord_y(jj) ) .AND. ( j < nyng ) )
              j  = j + 1
           ENDDO
           jflo(jj) = MIN( MAX( j, nysg ), nyng )
           jr = j
-          DO  WHILE ( ( coord_y(jr) + 0.5_wp * dy <= cg%coord_y(jj) + cg%dy )  &
-                      .AND.  ( j < nyng+1 ) )
+          DO WHILE ( ( coord_y(jr) + 0.5_wp * dy <= pg%coord_y(jj) + pg%dy ) .AND. ( j < nyng+1 ) )
              j  = j + 1
              jr = MIN( j, nyng )
@@ -1814,9 +1880,9 @@
 !
 !--       Print out the index bounds for checking and debugging purposes
-          WRITE(9,"('pmci_define_index_mapping, jj, jflo, jfuo: ', 3(i4,2x))")    &
+          WRITE( 9, "('pmci_define_index_mapping, jj, jflo, jfuo: ', 3(i4,2x))" )                   &
                jj, jflo(jj), jfuo(jj)
-          FLUSH(9)
+          FLUSH( 9 )
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    k-indices of w for each kk-index value
@@ -1826,10 +1892,13 @@
        kflw(0) = 0
        kfuw(0) = 0
-       DO kk = 1, cg%nz+1
+       DO  kk = 1, pg%nz+1
 !
 !--       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.
+!--       However, kctw still has to be stored separately as these index bounds 
+!--       are passed as arguments to the interpolation and anterpolation 
+!--       subroutines.
           k = kstart
-          DO WHILE ( ( zw(k) < cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
+          DO WHILE ( ( zw(k) < pg%zw(kk) )  .AND.  ( k < nzt+1 ) )
              k = k + 1
           ENDDO
@@ -1839,9 +1908,9 @@
 !
 !--       Print out the index bounds for checking and debugging purposes
-          WRITE(9,"('pmci_define_index_mapping, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))") &
-               kk, kflw(kk), kfuw(kk), nzt,  cg%zu(kk), cg%zw(kk)
-          FLUSH(9)
+          WRITE( 9, "('pmci_define_index_mapping, kk, kflw, kfuw: ', 4(i4,2x), 2(e12.5,2x))" )      &
+               kk, kflw(kk), kfuw(kk), nzt,  pg%zu(kk), pg%zw(kk)
+          FLUSH( 9 )
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    k-indices of others for each kk-index value
@@ -1852,11 +1921,11 @@
 !--    Note that anterpolation index limits are needed also for the top boundary 
 !--    ghost cell level because they are used also in the interpolation.
-       DO  kk = 1, cg%nz+1
+       DO  kk = 1, pg%nz+1
           k = kstart
-          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k <= nzt ) )
+          DO WHILE ( ( zu(k) < pg%zw(kk-1) )  .AND.  ( k <= nzt ) )
              k = k + 1
           ENDDO
           kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
+          DO WHILE ( ( zu(k) <= pg%zw(kk) )  .AND.  ( k <= nzt+1 ) )
              k = k + 1
              IF ( k > nzt + 1 ) EXIT  ! This EXIT is to prevent zu(k) from flowing over.
@@ -1866,26 +1935,28 @@
        ENDDO
 !
-!--    Set the k-index bounds separately for the parent-grid cells cg%nz and cg%nz+1       
+!--    Set the k-index bounds separately for the parent-grid cells pg%nz and pg%nz+1       
 !--    although they are not actually needed.
-       kflo(cg%nz)   = nzt+1    
-       kfuo(cg%nz)   = nzt+kgsr 
-       kflo(cg%nz+1) = nzt+kgsr 
-       kfuo(cg%nz+1) = nzt+kgsr 
+!  WHY IS THIS LIKE THIS? REVISE WITH CARE.
+       kflo(pg%nz)   = nzt+1    
+       kfuo(pg%nz)   = nzt+kgsr 
+       kflo(pg%nz+1) = nzt+kgsr 
+       kfuo(pg%nz+1) = nzt+kgsr 
 !
 !--    Print out the index bounds for checking and debugging purposes
-       DO  kk = 1, cg%nz+1
-          WRITE(9,"('pmci_define_index_mapping, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))") &
-               kk, kflo(kk), kfuo(kk), nzt,  cg%zu(kk), cg%zw(kk)
-          FLUSH(9)
+       DO  kk = 1, pg%nz+1
+          WRITE( 9, "('pmci_define_index_mapping, kk, kflo, kfuo: ', 4(i4,2x), 2(e12.5,2x))" )      &
+               kk, kflo(kk), kfuo(kk), nzt,  pg%zu(kk), pg%zw(kk)
+          FLUSH( 9 )
        ENDDO
-       WRITE(9,*)
+       WRITE( 9, * )
 !
 !--    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 and in reversibility 
-!--    correction in interpolation. 
-       DO  ii = icla, icra
-          DO  jj = jcsa, jcna
-             DO kk = 0, cg%nz+1
+!--    This information is needed in the anterpolation. 
+!--    The indices for wall_flags_0 (kw,jw,iw) must be limited to the range 
+!--    [-1,...,nx/ny/nzt+1] in order to avoid zero values on the outer ghost nodes.
+       DO  ii = ipla, ipra
+          DO  jj = jpsa, jpna
+             DO  kk = 0, pg%nz+1
 !
 !--             u-component
@@ -1896,5 +1967,5 @@
                       DO  k = kflo(kk), kfuo(kk)
                          kw = MIN( k, nzt+1 )                
-                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii)                  &
+                         ijkfc_u(kk,jj,ii) = ijkfc_u(kk,jj,ii)                                      &
                               + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 1 ) )
                       ENDDO
@@ -1909,5 +1980,5 @@
                       DO  k = kflo(kk), kfuo(kk)
                          kw = MIN( k, nzt+1 )                                        
-                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii)                  &
+                         ijkfc_v(kk,jj,ii) = ijkfc_v(kk,jj,ii)                                      &
                               + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 2 ) )
                       ENDDO
@@ -1922,5 +1993,5 @@
                       DO  k = kflo(kk), kfuo(kk)
                          kw = MIN( k, nzt+1 )
-                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii)                  &
+                         ijkfc_s(kk,jj,ii) = ijkfc_s(kk,jj,ii)                                      &
                               + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 0 ) )
                       ENDDO
@@ -1935,6 +2006,6 @@
                       DO  k = kflw(kk), kfuw(kk)
                          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 ) )
+                         ijkfc_w(kk,jj,ii) = ijkfc_w(kk,jj,ii)                                      &
+                              + MERGE( 1, 0, BTEST( wall_flags_0(kw,jw,iw), 3 ) )
                       ENDDO
                    ENDDO
@@ -1956,41 +2027,36 @@
 !
 !--    Determine and store the PE-subdomain dependent index bounds
-       IF  ( bc_dirichlet_l )  THEN
-          iclw = icl + 1
+       IF ( bc_dirichlet_l )  THEN
+          iplw = ipl + 1
        ELSE
-          iclw = icl - 1
-       ENDIF
-
-       IF  ( bc_dirichlet_r )  THEN
-          icrw = icr - 1
+          iplw = ipl - 1
+       ENDIF
+
+       IF ( bc_dirichlet_r )  THEN
+          iprw = ipr - 1
        ELSE
-          icrw = icr + 1
-       ENDIF
-
-       IF  ( bc_dirichlet_s )  THEN
-          jcsw = jcs + 1
+          iprw = ipr + 1
+       ENDIF
+
+       IF ( bc_dirichlet_s )  THEN
+          jpsw = jps + 1
        ELSE
-          jcsw = jcs - 1
-       ENDIF
-
-       IF  ( bc_dirichlet_n )  THEN
-          jcnw = jcn - 1
+          jpsw = jps - 1
+       ENDIF
+
+       IF ( bc_dirichlet_n )  THEN
+          jpnw = jpn - 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
+          jpnw = jpn + 1
+       ENDIF
 !
 !--    Left and right boundaries.
-       ALLOCATE( workarrc_lr(0:cg%nz+1,jcsw:jcnw,0:2) )
+       ALLOCATE( workarr_lr(0:pg%nz+1,jpsw:jpnw,0:2) )
 !
 !--    South and north boundaries.
-       ALLOCATE( workarrc_sn(0:cg%nz+1,0:2,iclw:icrw) )
+       ALLOCATE( workarr_sn(0:pg%nz+1,0:2,iplw:iprw) )
 !
 !--    Top boundary.
-       ALLOCATE( workarrc_t(0:2,jcsw:jcnw,iclw:icrw) )
+       ALLOCATE( workarr_t(0:2,jpsw:jpnw,iplw:iprw) )
 
     END SUBROUTINE pmci_allocate_workarrays
@@ -2000,29 +2066,27 @@
     SUBROUTINE pmci_create_workarray_exchange_datatypes
 !
-!--    Define specific MPI types for workarrc-exhchange.
+!--    Define specific MPI types for workarr-exchange.
        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 )
+       CALL MPI_TYPE_VECTOR( 3, pg%nz+2, (jpnw-jpsw+1)*(pg%nz+2), MPI_REAL,                         &
+            workarr_lr_exchange_type, ierr )
+       CALL MPI_TYPE_COMMIT( workarr_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 )
+       CALL MPI_TYPE_VECTOR( 1, 3*(pg%nz+2), 3*(pg%nz+2), MPI_REAL,                                 &
+            workarr_sn_exchange_type, ierr )
+       CALL MPI_TYPE_COMMIT( workarr_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 )
+       CALL MPI_TYPE_VECTOR( iprw-iplw+1, 3, 3*(jpnw-jpsw+1), MPI_REAL,                             &
+            workarr_t_exchange_type_x, ierr )
+       CALL MPI_TYPE_COMMIT( workarr_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
+       CALL MPI_TYPE_VECTOR( 1, 3*(jpnw-jpsw+1), 3*(jpnw-jpsw+1), MPI_REAL,                         &
+            workarr_t_exchange_type_y, ierr )
+       CALL MPI_TYPE_COMMIT( workarr_t_exchange_type_y, ierr )
        
     END SUBROUTINE pmci_create_workarray_exchange_datatypes
@@ -2036,96 +2100,115 @@
 !--    the parent grid spacing in the respective direction. 
        IMPLICIT NONE
-       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
-
-
-       non_matching_lower_left_corner = 0
-       non_int_gsr_x = 0
-       non_int_gsr_y = 0
-       non_int_gsr_z = 0
-       too_narrow_pesd_x = 0
-       too_narrow_pesd_y = 0
-
-       IF  ( myid == 0 )  THEN
+
+       INTEGER(iwp) ::  non_matching_height = 0              !< Flag for non-matching child-domain height
+       INTEGER(iwp) ::  non_matching_lower_left_corner = 0   !< Flag for non-matching lower left corner
+       INTEGER(iwp) ::  non_matching_upper_right_corner = 0  !< Flag for non-matching upper right corner
+       INTEGER(iwp) ::  non_int_gsr_x = 0                    !< Flag for non-integer grid-spacing ration in x-direction
+       INTEGER(iwp) ::  non_int_gsr_y = 0                    !< Flag for non-integer grid-spacing ration in y-direction
+       INTEGER(iwp) ::  non_int_gsr_z = 0                    !< Flag for non-integer grid-spacing ration in z-direction
+       INTEGER(iwp) ::  too_narrow_pesd_x = 0                !< Flag for too narrow pe-subdomain in x-direction
+       INTEGER(iwp) ::  too_narrow_pesd_y = 0                !< Flag for too narrow pe-subdomain in y-direction
+                                                         
+       REAL(wp), PARAMETER ::  tolefac = 1.0E-6_wp           !< Relative tolerence for grid-line matching
+                                                         
+       REAL(wp) ::  child_ngp_x_l                            !< Number of gridpoints in child subdomain in x-direction
+                                                             !< converted to REAL(wp) 
+       REAL(wp) ::  child_ngp_y_l                            !< Number of gridpoints in child subdomain in y-direction
+                                                             !< converted to REAL(wp)
+       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
+       REAL(wp) ::  upper_right_coord_x                      !< X-coordinate of the upper right corner of the child domain
+       REAL(wp) ::  upper_right_coord_y                      !< Y-coordinate of the upper right corner of the child domain
+
+       
+       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
+          tolez = tolefac * dz(1)
+!
+!--       First check that the child domain lower left corner matches the parent grid lines.
+          IF ( MOD( lower_left_coord_x, pg%dx ) > tolex ) non_matching_lower_left_corner = 1
+          IF ( MOD( lower_left_coord_y, pg%dy ) > toley ) non_matching_lower_left_corner = 1
+!
+!--       Then check that the child doman upper right corner matches the parent grid lines.
+          upper_right_coord_x = lower_left_coord_x + ( nx + 1 ) * dx
+          upper_right_coord_y = lower_left_coord_y + ( ny + 1 ) * dy
+          IF ( MOD( upper_right_coord_x, pg%dx ) > tolex ) non_matching_upper_right_corner = 1
+          IF ( MOD( upper_right_coord_y, pg%dy ) > toley ) non_matching_upper_right_corner = 1
+!
+!--       Also check that the cild domain height matches the parent grid lines.
+          IF ( MOD( zw(nzt), pg%dz ) > tolez ) non_matching_height = 1
+!
+!--       Check that the grid-spacing ratios in each direction are integer valued.    
+          IF ( MOD( pg%dx, dx ) > tolex )  non_int_gsr_x = 1
+          IF ( MOD( pg%dy, dy ) > toley )  non_int_gsr_y = 1
+!
+!--       In the z-direction, all levels need to be checked separately against grid stretching  
+!--       which is not allowed.
           DO  n = 0, kctw+1
-             IF  ( ABS( cg%zw(n) - zw(kflw(n)) ) > tolez )  non_int_gsr_z = 1
+             IF ( ABS( pg%zw(n) - zw(kflw(n)) ) > tolez )  non_int_gsr_z = 1
           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
-
-       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    
-
-       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
+          child_ngp_x_l = REAL( nxr - nxl + 1, KIND=wp )
+          IF ( child_ngp_x_l / REAL( igsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_x = 1
+          child_ngp_y_l = REAL( nyn - nys + 1, KIND=wp )
+          IF ( child_ngp_y_l / REAL( jgsr, KIND=wp ) < 1.0_wp )  too_narrow_pesd_y = 1
+         
+          IF ( non_matching_height > 0 )  THEN
+             WRITE( message_string, * ) 'nested child domain height must match ',                   &
+                                        'its parent grid lines'
+             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
+          ENDIF
+
+          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
+
+          IF ( non_matching_upper_right_corner > 0 )  THEN
+             WRITE( message_string, * ) 'nested child domain upper right ',                         &
+                                        'corner must match its parent grid lines'
+             CALL message( 'pmci_check_grid_matching', 'PA0414', 3, 2, 0, 6, 0 )
+          ENDIF
+
+          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
+
+          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
+
+          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
+
+          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
-
+          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
+                 
+       ENDIF  !  ( myid == 0 )
+       
+    END SUBROUTINE pmci_check_grid_matching
+     
 #endif  
  END SUBROUTINE pmci_setup_child
@@ -2203,14 +2286,11 @@
 !
 !-- Chemistry, depends on number of species
-    IF ( air_chemistry  .AND.  nest_chemistry )                                &
-       pmc_max_array = pmc_max_array + nspec
-!
-!-- Salsa, depens on the number aerosol size bins and chemical components + 
+    IF ( air_chemistry  .AND.  nest_chemistry )  pmc_max_array = pmc_max_array + nspec
+!
+!-- SALSA, depens on the number aerosol size bins and chemical components + 
 !-- the number of default gases
-    IF ( salsa  .AND.  nest_salsa )                                            &
-       pmc_max_array = pmc_max_array + nbins_aerosol + nbins_aerosol *         &
-                       ncomponents_mass
-       IF ( .NOT. salsa_gases_from_chem )  pmc_max_array = pmc_max_array +     &
-                                                           ngases_salsa
+    IF ( salsa  .AND.  nest_salsa )  pmc_max_array = pmc_max_array + nbins_aerosol +                &
+         nbins_aerosol * ncomponents_mass
+    IF ( .NOT. salsa_gases_from_chem )  pmc_max_array = pmc_max_array + ngases_salsa
 
 
@@ -2220,24 +2300,24 @@
 
 
-
- SUBROUTINE pmci_set_array_pointer( name, child_id, nz_cl, n )
+ SUBROUTINE pmci_set_array_pointer( name, child_id, nz_child, n )
 
     IMPLICIT NONE
-
-    INTEGER(iwp), INTENT(IN)          ::  child_id    !<
-    INTEGER(iwp), INTENT(IN)          ::  nz_cl       !<
-    INTEGER(iwp), INTENT(IN),OPTIONAL ::  n           !< index of chemical 
-                                                      !< species / salsa variables
-
+   
+    INTEGER(iwp), INTENT(IN) ::  child_id  !<
+    INTEGER(iwp), INTENT(IN) ::  nz_child  !<
+    
+    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n          !< index of chemical species
+    
     CHARACTER(LEN=*), INTENT(IN) ::  name             !<
-
-#if defined( __parallel )
-    INTEGER(iwp) ::  ierr                            !< MPI error code
-
-    REAL(wp), POINTER, DIMENSION(:,:)     ::  p_2d        !<
-    REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d        !<
-    REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d_sec    !<
-    INTEGER(idp), POINTER, DIMENSION(:,:) ::  i_2d        !<
-
+!
+!-- Local variables:       
+    INTEGER(iwp) ::  ierr                             !< MPI error code
+
+    INTEGER(idp), POINTER, DIMENSION(:,:) ::  i_2d    !<
+       
+    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d      !<
+    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d      !<
+    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d_sec  !<
+    
 
     NULLIFY( p_3d )
@@ -2267,6 +2347,6 @@
        p_3d => salsa_gas(n)%conc
 !
-!-- Next line is just an example for a 2D array (not active for coupling!)
-!-- Please note, that z0 has to be declared as TARGET array in modules.f90
+!-- Next line is just an example for a 2D array (not active for coupling!) 
+!-- Please note, that z0 has to be declared as TARGET array in modules.f90.
 !    IF ( TRIM(name) == "z0" )    p_2d => z0
     IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
@@ -2289,17 +2369,14 @@
 
     IF ( ASSOCIATED( p_3d ) )  THEN
-       CALL pmc_s_set_dataarray( child_id, p_3d, nz_cl, nz,                    &
-                                 array_2 = p_3d_sec )
-    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
+       CALL pmc_s_set_dataarray( child_id, p_3d, nz_child, nz, array_2 = p_3d_sec )
+    ELSEIF  ( ASSOCIATED( p_2d ) )  THEN
        CALL pmc_s_set_dataarray( child_id, p_2d )
-    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
+    ELSEIF  ( ASSOCIATED( i_2d ) )  THEN
        CALL pmc_s_set_dataarray( child_id, i_2d )
     ELSE
 !
 !--    Give only one message for the root domain
-       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
-
-          message_string = 'pointer for array "' // TRIM( name ) //            &
-                           '" can''t be associated'
+       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN             
+          message_string = 'pointer for array "' // TRIM( name ) // '" can''t be associated'
           CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
        ELSE
@@ -2308,112 +2385,120 @@
           CALL MPI_BARRIER( comm2d, ierr )
        ENDIF
-
-   ENDIF
-
+       
+    ENDIF
+    
+ END SUBROUTINE pmci_set_array_pointer
+
+
+     
+ INTEGER FUNCTION get_number_of_children()
+
+    IMPLICIT NONE
+
+    
+#if defined( __parallel )
+    get_number_of_children = SIZE( pmc_parent_for_child ) - 1
+#else
+    get_number_of_children = 0
 #endif
-END SUBROUTINE pmci_set_array_pointer
-
-
-
-INTEGER FUNCTION get_number_of_childs ()     ! Change the name to "get_number_of_children"
-
-   IMPLICIT NONE
-
+
+    RETURN
+
+ END FUNCTION get_number_of_children
+
+
+ 
+ INTEGER FUNCTION get_childid( id_index )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) ::  id_index   !<
+
+    
 #if defined( __parallel )
-   get_number_of_childs = SIZE( pmc_parent_for_child ) - 1
+    get_childid = pmc_parent_for_child(id_index)
 #else
-   get_number_of_childs = 0
+    get_childid = 0
 #endif
 
-   RETURN
-
-END FUNCTION get_number_of_childs
-
-
-INTEGER FUNCTION get_childid (id_index)
-
-   IMPLICIT NONE
-
-   INTEGER,INTENT(IN)                 :: id_index
-
-#if defined( __parallel )
-   get_childid = pmc_parent_for_child(id_index)
-#else
-   get_childid = 0
-#endif
-
-   RETURN
-
-END FUNCTION get_childid
-
-
-
-SUBROUTINE  get_child_edges (m, lx_coord, lx_coord_b, rx_coord, rx_coord_b,    &
-                               sy_coord, sy_coord_b, ny_coord, ny_coord_b,     &
-                               uz_coord, uz_coord_b)
-   IMPLICIT NONE
-   INTEGER,INTENT(IN)             ::  m
-   REAL(wp),INTENT(OUT)           ::  lx_coord, lx_coord_b
-   REAL(wp),INTENT(OUT)           ::  rx_coord, rx_coord_b
-   REAL(wp),INTENT(OUT)           ::  sy_coord, sy_coord_b
-   REAL(wp),INTENT(OUT)           ::  ny_coord, ny_coord_b
-   REAL(wp),INTENT(OUT)           ::  uz_coord, uz_coord_b
-
-   lx_coord = childgrid(m)%lx_coord
-   rx_coord = childgrid(m)%rx_coord
-   sy_coord = childgrid(m)%sy_coord
-   ny_coord = childgrid(m)%ny_coord
-   uz_coord = childgrid(m)%uz_coord
-
-   lx_coord_b = childgrid(m)%lx_coord_b
-   rx_coord_b = childgrid(m)%rx_coord_b
-   sy_coord_b = childgrid(m)%sy_coord_b
-   ny_coord_b = childgrid(m)%ny_coord_b
-   uz_coord_b = childgrid(m)%uz_coord_b
-
-END SUBROUTINE get_child_edges
-
-
-
-SUBROUTINE  get_child_gridspacing( m, dx,dy,dz )
-
-   IMPLICIT NONE
-   INTEGER,INTENT(IN)             ::  m
-   REAL(wp),INTENT(OUT)           ::  dx,dy
-   REAL(wp),INTENT(OUT),OPTIONAL  ::  dz
-
-   dx = childgrid(m)%dx
-   dy = childgrid(m)%dy
-   IF(PRESENT(dz))   THEN
-      dz = childgrid(m)%dz
-   ENDIF
-
-END SUBROUTINE get_child_gridspacing
-
-
-
-SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc, n  )
-
+    RETURN
+
+ END FUNCTION get_childid
+
+
+
+ SUBROUTINE get_child_edges( m, lx_coord, lx_coord_b, rx_coord, rx_coord_b, sy_coord, sy_coord_b,   &
+      ny_coord, ny_coord_b, uz_coord, uz_coord_b )
+   
     IMPLICIT NONE
 
-    CHARACTER(LEN=*), INTENT(IN) ::  name    !<
-
-    INTEGER(iwp), INTENT(IN) ::  ie      !<
+    INTEGER,INTENT(IN)   ::  m                     !<
+
+    REAL(wp),INTENT(OUT) ::  lx_coord, lx_coord_b  !<
+    REAL(wp),INTENT(OUT) ::  rx_coord, rx_coord_b  !<
+    REAL(wp),INTENT(OUT) ::  sy_coord, sy_coord_b  !<
+    REAL(wp),INTENT(OUT) ::  ny_coord, ny_coord_b  !<
+    REAL(wp),INTENT(OUT) ::  uz_coord, uz_coord_b  !<
+
+    
+    lx_coord = childgrid(m)%lx_coord
+    rx_coord = childgrid(m)%rx_coord
+    sy_coord = childgrid(m)%sy_coord
+    ny_coord = childgrid(m)%ny_coord
+    uz_coord = childgrid(m)%uz_coord
+    
+    lx_coord_b = childgrid(m)%lx_coord_b
+    rx_coord_b = childgrid(m)%rx_coord_b
+    sy_coord_b = childgrid(m)%sy_coord_b
+    ny_coord_b = childgrid(m)%ny_coord_b
+    uz_coord_b = childgrid(m)%uz_coord_b
+    
+ END SUBROUTINE get_child_edges
+
+
+
+ SUBROUTINE  get_child_gridspacing( m, dx, dy,dz )
+
+    IMPLICIT NONE
+   
+    INTEGER, INTENT(IN)             ::  m      !<
+
+    REAL(wp), INTENT(OUT)           ::  dx,dy  !<
+
+    REAL(wp), INTENT(OUT), OPTIONAL ::  dz     !<
+
+    
+    dx = childgrid(m)%dx
+    dy = childgrid(m)%dy
+    IF ( PRESENT( dz ) )  THEN
+       dz = childgrid(m)%dz
+    ENDIF
+    
+ END SUBROUTINE get_child_gridspacing
+
+
+
+ SUBROUTINE pmci_create_childs_parent_grid_arrays( name, is, ie, js, je, nzc, n  )
+
+    IMPLICIT NONE
+
+    INTEGER(iwp), INTENT(IN) ::  ie      !<  RENAME ie, is, je, js?
     INTEGER(iwp), INTENT(IN) ::  is      !<
     INTEGER(iwp), INTENT(IN) ::  je      !<
     INTEGER(iwp), INTENT(IN) ::  js      !<
-    INTEGER(iwp), INTENT(IN) ::  nzc     !<  nzc is cg%nz, but note that cg%nz is not the original nz of parent, but the highest parent-grid level needed for nesting.
-
-    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species /
-                                              !< salsa variables
-
-#if defined( __parallel )
+    INTEGER(iwp), INTENT(IN) ::  nzc     !<  nzc is pg%nz, but note that pg%nz is not the original nz of parent,
+                                            !<  but the highest parent-grid level needed for nesting.
+    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species / salsa variables
+    
+    CHARACTER(LEN=*), INTENT(IN) ::  name    !<
+!       
+!-- Local variables:
     INTEGER(iwp) ::  ierr    !<
-
+    
+    INTEGER(idp), POINTER,DIMENSION(:,:)   ::  i_2d    !<
+    
     REAL(wp), POINTER,DIMENSION(:,:)       ::  p_2d    !<
     REAL(wp), POINTER,DIMENSION(:,:,:)     ::  p_3d    !<
-    INTEGER(idp), POINTER,DIMENSION(:,:)   ::  i_2d    !<
-
-
+    
     NULLIFY( p_3d )
     NULLIFY( p_2d )
@@ -2464,6 +2549,5 @@
        i_2d => part_adrc
     ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
-       IF ( .NOT. ALLOCATED( chem_spec_c ) )                                   &
-          ALLOCATE( chem_spec_c(0:nzc+1,js:je,is:ie,1:nspec) )
+       IF ( .NOT. ALLOCATED( chem_spec_c ) ) ALLOCATE( chem_spec_c(0:nzc+1,js:je,is:ie,1:nspec) )
        p_3d => chem_spec_c(:,:,:,n)
     ELSEIF ( TRIM( name(1:3) ) == "an_" )  THEN
@@ -2493,22 +2577,23 @@
     ELSE
 !
-!--    Give only one message for the first child domain
+!-- Give only one message for the first child domain.
        IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
-
           message_string = 'pointer for array "' // TRIM( name ) //            &
-                           '" can''t be associated'
-          CALL message( 'pmci_create_child_arrays', 'PA0170', 3, 2, 0, 6, 0 )
+               '" can''t be associated'
+          CALL message( 'pmci_create_childs_parent_grid_arrays', 'PA0170', 3, 2, 0, 6, 0 )
        ELSE
 !
-!--       Prevent others from continuing
+!--          Prevent others from continuing in case the abort is to come.
           CALL MPI_BARRIER( comm2d, ierr )
        ENDIF
+
     ENDIF
 
-#endif
- END SUBROUTINE pmci_create_child_arrays
-
-
-
+ END SUBROUTINE pmci_create_childs_parent_grid_arrays
+
+
+!
+! E N D   O F    S E T U P   R O U T I N E S 
+!
  SUBROUTINE pmci_parent_initialize
 
@@ -2542,44 +2627,16 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i          !<
-    INTEGER(iwp) ::  ib         !< running index for aerosol size bins
-    INTEGER(iwp) ::  ic         !< running index for aerosol mass bins
-    INTEGER(iwp) ::  icl        !<
-    INTEGER(iwp) ::  icla       !<
-    INTEGER(iwp) ::  iclw       !<
-    INTEGER(iwp) ::  icr        !<
-    INTEGER(iwp) ::  icra       !<
-    INTEGER(iwp) ::  icrw       !<
-    INTEGER(iwp) ::  ig         !< running index for salsa gases
-    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
-
-    REAL(wp) ::  waittime       !<
+    INTEGER(iwp) ::  ic         !< Child-grid index in x-direction
+    INTEGER(iwp) ::  jc         !< Child-grid index in y-direction
+    INTEGER(iwp) ::  kc         !< Child-grid index in z-direction
+    INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
+    INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins 
+    INTEGER(iwp) ::  lg         !< Running index for salsa gases
+    INTEGER(iwp) ::  n          !< Running index for chemical species
+    REAL(wp) ::  waittime       !< Waiting time 
 
 !
 !-- Root model is never anyone's child
     IF ( cpl_id > 1 )  THEN
-!
-!--    Child domain boundaries in the parent index space
-       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
@@ -2633,15 +2690,15 @@
 
        IF ( salsa  .AND.  nest_salsa )  THEN
-          DO  ib = 1, nbins_aerosol
-             CALL pmci_interp_1sto_all ( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),        &
+          DO  lb = 1, nbins_aerosol
+             CALL pmci_interp_1sto_all ( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),       &
                                          kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
           ENDDO
-          DO  ic = 1, nbins_aerosol * ncomponents_mass
-             CALL pmci_interp_1sto_all ( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),            &
+          DO  lc = 1, nbins_aerosol * ncomponents_mass
+             CALL pmci_interp_1sto_all ( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),           &
                                          kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
           ENDDO
           IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                CALL pmci_interp_1sto_all ( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),             &
+             DO  lg = 1, ngases_salsa
+                CALL pmci_interp_1sto_all ( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),              &
                                             kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 's' )
              ENDDO
@@ -2651,26 +2708,14 @@
        IF ( topography /= 'flat' )  THEN
 !
-!--       Inside buildings set velocities and TKE back to zero.
-!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
-!--       maybe revise later.
-          DO   i = nxlg, nxrg
-             DO   j = nysg, nyng
-                DO  k = nzb, nzt
-                   u(k,j,i)   = MERGE( u(k,j,i), 0.0_wp,                       &
-                                       BTEST( wall_flags_0(k,j,i), 1 ) )
-                   v(k,j,i)   = MERGE( v(k,j,i), 0.0_wp,                       &
-                                       BTEST( wall_flags_0(k,j,i), 2 ) )
-                   w(k,j,i)   = MERGE( w(k,j,i), 0.0_wp,                       &
-                                       BTEST( wall_flags_0(k,j,i), 3 ) )
-!                    e(k,j,i)   = MERGE( e(k,j,i), 0.0_wp,                       &
-!                                        BTEST( wall_flags_0(k,j,i), 0 ) )
-                   u_p(k,j,i) = MERGE( u_p(k,j,i), 0.0_wp,                     &
-                                       BTEST( wall_flags_0(k,j,i), 1 ) )
-                   v_p(k,j,i) = MERGE( v_p(k,j,i), 0.0_wp,                     &
-                                       BTEST( wall_flags_0(k,j,i), 2 ) )
-                   w_p(k,j,i) = MERGE( w_p(k,j,i), 0.0_wp,                     &
-                                       BTEST( wall_flags_0(k,j,i), 3 ) )
-!                    e_p(k,j,i) = MERGE( e_p(k,j,i), 0.0_wp,                     &
-!                                        BTEST( wall_flags_0(k,j,i), 0 ) )
+!--       Inside buildings set velocities back to zero.
+          DO  ic = nxlg, nxrg
+             DO  jc = nysg, nyng
+                DO  kc = nzb, nzt
+                   u(kc,jc,ic)   = MERGE( u(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 1 ) )
+                   v(kc,jc,ic)   = MERGE( v(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 2 ) )
+                   w(kc,jc,ic)   = MERGE( w(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 3 ) )
+                   u_p(kc,jc,ic) = MERGE( u_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 1 ) )
+                   v_p(kc,jc,ic) = MERGE( v_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 2 ) )
+                   w_p(kc,jc,ic) = MERGE( w_p(kc,jc,ic), 0.0_wp, BTEST( wall_flags_0(kc,jc,ic), 3 ) )
                 ENDDO
              ENDDO
@@ -2683,97 +2728,106 @@
 
 
-    SUBROUTINE pmci_interp_1sto_all( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu, var )
+    SUBROUTINE pmci_interp_1sto_all( child_array, parent_array, kct, ifl, ifu, jfl, jfu, kfl, kfu,  &
+         var )
 !
 !--    Interpolation of the internal values for the child-domain initialization
        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
-       INTEGER(iwp) :: kct                                    !< The parent-grid index in z-direction just below the boundary value node
-       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 end 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 end index of child cells belonging to certain parent cell - z direction
-       CHARACTER(LEN=1), INTENT(IN) :: var                    !< Variable symbol: 'u', 'v', 'w' or 's'
+       INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
+
+       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !<  Indicates start index of child cells belonging to certain
+                                                               !<  parent cell - x direction
+       INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !<  Indicates end index of child cells belonging to certain
+                                                               !<  parent cell - x direction
+       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !<  Indicates start index of child cells belonging to certain
+                                                               !<  parent cell - y direction
+       INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !<  Indicates end index of child cells belonging to certain
+                                                               !<  parent cell - y direction
+       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !<  Indicates start index of child cells belonging to certain
+                                                               !<  parent cell - z direction
+       INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !<  Indicates end index of child cells belonging to certain
+                                                               !<  parent cell - z direction
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array  !<  Child-grid array 
+       REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN) ::  parent_array        !<  Parent-grid array
+
+       CHARACTER(LEN=1), INTENT(IN) ::  var  !<  Variable symbol: 'u', 'v', 'w' or 's'
 !
 !--    Local variables:
-       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        !<
-       INTEGER(iwp) ::  lb       !<
-       INTEGER(iwp) ::  le       !<
-       INTEGER(iwp) ::  m        !<
-       INTEGER(iwp) ::  mb       !<
-       INTEGER(iwp) ::  me       !<
-       INTEGER(iwp) ::  n        !<
-
-
-       lb = icl
-       le = icr
-       mb = jcs
-       me = jcn
-       ifirst = nxl
-       ilast  = nxr
-       jfirst = nys
-       jlast  = nyn
+       INTEGER(iwp) ::  ic        !< Running child-grid index in the x-direction
+       INTEGER(iwp) ::  icfirst   !< Leftmost child-grid index initialized by the main loops (usually icfirst == icl_init)
+       INTEGER(iwp) ::  iclast    !< Rightmost child-grid index initialized by the main loops (usually iclast == icr_init)
+       INTEGER(iwp) ::  icl_init  !< Left child-grid index bound for initialization in the x-direction
+       INTEGER(iwp) ::  icr_init  !< Right child-grid index bound for initialization in the x-direction
+       INTEGER(iwp) ::  jc        !< Running child-grid index in the y-direction
+       INTEGER(iwp) ::  jcfirst   !< Southmost child-grid index initialized by the main loops (usually jcfirst == jcs_init)
+       INTEGER(iwp) ::  jclast    !< Northmost child-grid index initialized by the main loops (usually jclast == jcn_init)
+       INTEGER(iwp) ::  jcs_init  !< South child-grid index bound for initialization in the y-direction
+       INTEGER(iwp) ::  jcn_init  !< North child-grid index bound for initialization in the y-direction
+       INTEGER(iwp) ::  kc        !< Running child-grid index in the z-direction
+       INTEGER(iwp) ::  ip        !< Running parent-grid index in the x-direction
+       INTEGER(iwp) ::  ipl_init  !< Left parent-grid index bound for initialization in the x-direction
+       INTEGER(iwp) ::  ipr_init  !< Right parent-grid index bound for initialization in the x-direction
+       INTEGER(iwp) ::  jp        !< Running parent-grid index in the y-direction
+       INTEGER(iwp) ::  jps_init  !< South parent-grid index bound for initialization in the y-direction
+       INTEGER(iwp) ::  jpn_init  !< North parent-grid index bound for initialization in the y-direction
+       INTEGER(iwp) ::  kp        !< Running parent-grid index in the z-direction
+
+
+       ipl_init = ipl
+       ipr_init = ipr
+       jps_init = jps
+       jpn_init = jpn
+       icl_init = nxl
+       icr_init = nxr
+       jcs_init = nys
+       jcn_init = nyn
 
        IF ( nesting_mode /= 'vertical' )  THEN
           IF ( bc_dirichlet_l )  THEN
-             lb     = icl + 1
-             ifirst = nxl - 1
+             ipl_init = ipl + 1
+             icl_init = nxl - 1
 !
 !--          For u, nxl is a ghost node, but not for the other variables
              IF ( var == 'u' )  THEN
-                lb     = icl + 2
-                ifirst = nxl
+                ipl_init = ipl + 2
+                icl_init = nxl
              ENDIF
           ENDIF
           IF ( bc_dirichlet_s )  THEN
-             mb     = jcs + 1
-             jfirst = nys - 1 
+             jps_init = jps + 1
+             jcs_init = nys - 1 
 !
 !--          For v, nys is a ghost node, but not for the other variables
              IF ( var == 'v' )  THEN
-                mb     = jcs + 2
-                jfirst = nys
+                jps_init = jps + 2
+                jcs_init = nys
              ENDIF
           ENDIF
           IF ( bc_dirichlet_r )  THEN
-             le    = icr - 1
-             ilast = nxr + 1 
+             ipr_init = ipr - 1
+             icr_init = nxr + 1 
           ENDIF
           IF ( bc_dirichlet_n )  THEN
-             me    = jcn - 1
-             jlast = nyn + 1
+             jpn_init = jpn - 1
+             jcn_init = nyn + 1
           ENDIF
        ENDIF      
 
-       f(:,:,:) = 0.0_wp
-
-       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
-                DO n = 0, kct + 1 
-
-                   DO  i = ifl(l), ifl(l+1)-1
-                      DO  j = jfl(m), jfu(m)
-                         DO  k = kfl(n), MIN( kfu(n), nzt+1 )
-                            f(k,j,i) = fc(n,m,l)
+       child_array(:,:,:) = 0.0_wp
+
+       IF ( var == 'u' )  THEN
+
+          icfirst = ifl(ipl_init) 
+          iclast  = ifl(ipr_init+1) - 1
+          jcfirst = jfl(jps_init)
+          jclast  = jfu(jpn_init)
+          DO  ip = ipl_init, ipr_init
+             DO  jp = jps_init, jpn_init
+                DO  kp = 0, kct + 1 
+
+                   DO  ic = ifl(ip), ifl(ip+1)-1
+                      DO  jc = jfl(jp), jfu(jp)
+                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
+                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
                          ENDDO
                       ENDDO
@@ -2784,18 +2838,18 @@
           ENDDO
 
-       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
-                DO n = 0, kct + 1  
-
-                   DO i = ifl(l), ifu(l)
-                      DO  j = jfl(m), jfl(m+1)-1
-                         DO  k = kfl(n), MIN( kfu(n), nzt+1 )
-                            f(k,j,i) = fc(n,m,l) 
+       ELSE IF ( var == 'v' )  THEN
+
+          icfirst = ifl(ipl_init) 
+          iclast  = ifu(ipr_init)
+          jcfirst = jfl(jps_init)
+          jclast  = jfl(jpn_init+1) - 1
+          DO  ip = ipl_init, ipr_init
+             DO  jp = jps_init, jpn_init
+                DO  kp = 0, kct + 1  
+
+                   DO  ic = ifl(ip), ifu(ip)
+                      DO  jc = jfl(jp), jfl(jp+1)-1
+                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
+                            child_array(kc,jc,ic) = parent_array(kp,jp,ip) 
                          ENDDO
                       ENDDO
@@ -2806,19 +2860,21 @@
           ENDDO
 
-       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
-                DO n = 1, 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)
+       ELSE IF ( var == 'w' )  THEN
+
+          icfirst = ifl(ipl_init) 
+          iclast  = ifu(ipr_init)
+          jcfirst = jfl(jps_init)
+          jclast  = jfu(jpn_init)
+          DO  ip = ipl_init, ipr_init
+             DO  jp = jps_init, jpn_init
+                DO  kp = 1, kct + 1  
+
+                   DO  ic = ifl(ip), ifu(ip)
+                      DO  jc = jfl(jp), jfu(jp)
+!                         
+!--                      Because the kp-loop for w starts from kp=1 instead of 0
+                         child_array(nzb,jc,ic) = 0.0_wp
+                         DO  kc = kfu(kp-1)+1, kfu(kp) 
+                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
                          ENDDO
                       ENDDO
@@ -2831,16 +2887,16 @@
        ELSE   ! scalars
 
-          ib = ifl(lb) 
-          ie = ifu(le)
-          jb = jfl(mb)
-          je = jfu(me)
-          DO  l = lb, le
-             DO  m = mb, me
-                DO n = 0, kct + 1
+          icfirst = ifl(ipl_init) 
+          iclast  = ifu(ipr_init)
+          jcfirst = jfl(jps_init)
+          jclast  = jfu(jpn_init)
+          DO  ip = ipl_init, ipr_init
+             DO  jp = jps_init, jpn_init
+                DO  kp = 0, kct + 1
                                       
-                   DO i = ifl(l), ifu(l)
-                      DO  j = jfl(m), jfu(m)                         
-                         DO  k = kfl(n), MIN( kfu(n), nzt+1 )
-                            f(k,j,i) = fc(n,m,l)
+                   DO  ic = ifl(ip), ifu(ip)
+                      DO  jc = jfl(jp), jfu(jp)                         
+                         DO  kc = kfl(kp), MIN( kfu(kp), nzt+1 )
+                            child_array(kc,jc,ic) = parent_array(kp,jp,ip)
                          ENDDO
                       ENDDO
@@ -2853,26 +2909,26 @@
        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)
+!--    If the number of grid points in child subdomain in x- or y-direction
+!--    (nxr - nxl + 1 and/or nyn - nys + 1) is not integer divisible by the grid spacing
+!--    ratio in its direction (igsr and/or jgsr), the above loops will return with 
+!--    unfilled gaps in the initial fields. These gaps, if present, are filled here.  
+       IF ( icfirst > icl_init )  THEN
+          DO  ic = icl_init, icfirst - 1
+             child_array(:,:,ic) = child_array(:,:,icfirst)
           ENDDO
        ENDIF
-       IF  ( ie < ilast )  THEN
-          DO  i = ie+1, ilast
-             f(:,:,i) = f(:,:,ie)
+       IF ( iclast < icr_init )  THEN
+          DO  ic = iclast + 1, icr_init
+             child_array(:,:,ic) = child_array(:,:,iclast)
           ENDDO
        ENDIF
-       IF  ( jb > jfirst )  THEN
-          DO  j = jfirst, jb-1
-             f(:,j,:) = f(:,jb,:)
+       IF ( jcfirst > jcs_init )  THEN
+          DO  jc = jcs_init, jcfirst - 1
+             child_array(:,jc,:) = child_array(:,jcfirst,:)
           ENDDO
        ENDIF
-       IF  ( je < jlast )  THEN
-          DO  j = je+1, jlast
-             f(:,j,:) = f(:,je,:)
+       IF ( jclast < jcn_init )  THEN
+          DO  jc = jclast + 1, jcn_init
+             child_array(:,jc,:) = child_array(:,jclast,:)
           ENDDO
        ENDIF
@@ -2887,5 +2943,5 @@
  SUBROUTINE pmci_check_setting_mismatches
 !
-!-- Check for mismatches between settings of master and child variables
+!-- Check for mismatches between settings of root and child variables
 !-- (e.g., all children have to follow the end_time settings of the root model).
 !-- The root model overwrites variables in the other models, so these variables
@@ -2894,15 +2950,15 @@
 #if defined( __parallel )
 
-    USE control_parameters,                                                    &
+    USE control_parameters,                                                                         &
         ONLY:  dt_restart, end_time, message_string, restart_time, time_restart
 
     IMPLICIT NONE
 
-    INTEGER ::  ierr
-
-    REAL(wp) ::  dt_restart_root
-    REAL(wp) ::  end_time_root
-    REAL(wp) ::  restart_time_root
-    REAL(wp) ::  time_restart_root
+    INTEGER ::  ierr                 !<  MPI error code
+
+    REAL(wp) ::  dt_restart_root     !< 
+    REAL(wp) ::  end_time_root       !<  
+    REAL(wp) ::  restart_time_root   !< 
+    REAL(wp) ::  time_restart_root   !<  
 
 !
@@ -2916,9 +2972,9 @@
     IF ( .NOT. pmc_is_rootmodel() )  THEN
        IF ( end_time /= end_time_root )  THEN
-          WRITE( message_string, * )  'mismatch between root model and ',      &
-               'child settings:& end_time(root) = ', end_time_root,            &
-               '& end_time(child) = ', end_time, '& child value is set',       &
+          WRITE( message_string, * )  'mismatch between root model and ',                           &
+               'child settings:& end_time(root) = ', end_time_root,                                 &
+               '& end_time(child) = ', end_time, '& child value is set',                            &
                ' to root value'
-          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
+          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6,                      &
                         0 )
           end_time = end_time_root
@@ -2990,12 +3046,12 @@
    IMPLICIT NONE
 
-   INTEGER(iwp)           :: ierr  !<
-   REAL(wp), DIMENSION(1) :: dtl   !<
-   REAL(wp), DIMENSION(1) :: dtg   !<
+   INTEGER(iwp) ::  ierr  !<  MPI error code
+   REAL(wp) ::  dtl       !<  Local time step of the current process
+   REAL(wp) ::  dtg       !<  Global time step defined as the global minimum of dtl of all processes
 
    
-   dtl(1) = dt_3d 
+   dtl = dt_3d
    CALL MPI_ALLREDUCE( dtl, dtg, 1, MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
-   dt_3d  = dtg(1)
+   dt_3d  = dtg
 
 #endif
@@ -3012,12 +3068,11 @@
     IMPLICIT NONE
 
-    INTEGER(iwp), INTENT(IN) ::  swaplevel  !< swaplevel (1 or 2) of PALM's
-                                            !< timestep
-
-    INTEGER(iwp)            ::  child_id    !<
-    INTEGER(iwp)            ::  m           !<
+    INTEGER(iwp), INTENT(IN) ::  swaplevel  !< swaplevel (1 or 2) of PALM's timestep
+
+    INTEGER(iwp) ::  child_id    !<  Child id of the child number m
+    INTEGER(iwp) ::  m           !<  Loop index over all children of the current parent 
 
 #if defined( __parallel )
-    DO  m = 1, SIZE( pmc_parent_for_child )-1
+    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
        child_id = pmc_parent_for_child(m)
        CALL pmc_s_set_active_data_array( child_id, swaplevel )
@@ -3042,5 +3097,5 @@
     IMPLICIT NONE
 
-    CHARACTER(LEN=*), INTENT(IN) ::  local_nesting_mode
+    CHARACTER(LEN=*), INTENT(IN) ::  local_nesting_mode  !<  Nesting mode: 'one-way', 'two-way' or 'vertical'
 
 
@@ -3052,5 +3107,5 @@
     ELSE
 
-       IF( nesting_datatransfer_mode == 'cascade' )  THEN
+       IF ( nesting_datatransfer_mode == 'cascade' )  THEN
 
           CALL pmci_child_datatrans( parent_to_child )
@@ -3060,5 +3115,5 @@
           CALL pmci_child_datatrans( child_to_parent )
 
-       ELSEIF( nesting_datatransfer_mode == 'overlap')  THEN
+       ELSEIF ( nesting_datatransfer_mode == 'overlap')  THEN
 
           CALL pmci_parent_datatrans( parent_to_child )
@@ -3068,5 +3123,5 @@
           CALL pmci_parent_datatrans( child_to_parent )
 
-       ELSEIF( TRIM( nesting_datatransfer_mode ) == 'mixed' )  THEN
+       ELSEIF ( TRIM( nesting_datatransfer_mode ) == 'mixed' )  THEN
 
           CALL pmci_parent_datatrans( parent_to_child )
@@ -3088,12 +3143,12 @@
     IMPLICIT NONE
 
-    INTEGER(iwp), INTENT(IN) ::  direction   !<
+    INTEGER(iwp), INTENT(IN) ::  direction   !<  Direction of the data transfer: 'parent_to_child' or 'child_to_parent'
 
 #if defined( __parallel )
-    INTEGER(iwp) ::  child_id    !<
-    INTEGER(iwp) ::  i           !<
-    INTEGER(iwp) ::  j           !<
-    INTEGER(iwp) ::  k           !<
-    INTEGER(iwp) ::  m           !<
+    INTEGER(iwp) ::  child_id    !<  Child id of the child number m
+    INTEGER(iwp) ::  i           !<  Parent-grid index in x-direction
+    INTEGER(iwp) ::  j           !<  Parent-grid index in y-direction
+    INTEGER(iwp) ::  k           !<  Parent-grid index in z-direction
+    INTEGER(iwp) ::  m           !<  Loop index over all children of the current parent 
 
 
@@ -3124,10 +3179,7 @@
                 DO   j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      u(k,j,i)  = MERGE( u(k,j,i), 0.0_wp,                     &
-                                         BTEST( wall_flags_0(k,j,i), 1 ) )
-                      v(k,j,i)  = MERGE( v(k,j,i), 0.0_wp,                     &
-                                         BTEST( wall_flags_0(k,j,i), 2 ) )
-                      w(k,j,i)  = MERGE( w(k,j,i), 0.0_wp,                     &
-                                         BTEST( wall_flags_0(k,j,i), 3 ) )
+                      u(k,j,i) = MERGE( u(k,j,i), 0.0_wp, BTEST( wall_flags_0(k,j,i), 1 ) )
+                      v(k,j,i) = MERGE( v(k,j,i), 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) )
+                      w(k,j,i) = MERGE( w(k,j,i), 0.0_wp, BTEST( wall_flags_0(k,j,i), 3 ) )
 !
 !--                 TO_DO: zero setting of temperature within topography creates
@@ -3142,5 +3194,5 @@
           ENDIF
        ENDIF
-    ENDDO
+    ENDDO  ! m
 
 #endif
@@ -3156,36 +3208,10 @@
 
 #if defined( __parallel )
-    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
+
+    REAL(wp), DIMENSION(1) ::  dtl          !< Time step size
 
 
     dtl = dt_3d
     IF ( cpl_id > 1 )  THEN
-!
-!--    Child domain boundaries in the parent indice space.
-       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
@@ -3213,5 +3239,5 @@
     ENDIF
 
-  CONTAINS
+ CONTAINS
 
    
@@ -3223,10 +3249,10 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  ibgp       !< index running over the nbgp boundary ghost points in i-direction
-       INTEGER(iwp) ::  jbgp       !< index running over the nbgp boundary ghost points in j-direction
-       INTEGER(iwp) ::  ib         !< running index for aerosol size bins
-       INTEGER(iwp) ::  ic         !< running index for aerosol mass bins
-       INTEGER(iwp) ::  ig         !< running index for salsa gases
-       INTEGER(iwp) ::  n          !< running index for number of chemical species
+       INTEGER(iwp) ::  ibgp       !< Index running over the nbgp boundary ghost points in i-direction
+       INTEGER(iwp) ::  jbgp       !< Index running over the nbgp boundary ghost points in j-direction
+       INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
+       INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
+       INTEGER(iwp) ::  lg         !< Running index for salsa gases
+       INTEGER(iwp) ::  n          !< Running index for number of chemical species
       
 !
@@ -3242,11 +3268,11 @@
              CALL pmci_interp_1sto_lr( w, wc, kctw, jflo, jfuo, kflw, kfuw, 'l', 'w' )
 
-             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                        &
-                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                         &
-                     .NOT. constant_diffusion ) )  THEN
+             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
+                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
+                     .NOT.  constant_diffusion ) )  THEN
 !                CALL pmci_interp_1sto_lr( e, ec, kcto, jflo, jfuo, kflo, kfuo, 'l', 'e' )
 !
 !--             Interpolation of e is replaced by the Neumann condition. 
-                DO ibgp = -nbgp, -1
+                DO  ibgp = -nbgp, -1
                    e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,0)
                 ENDDO
@@ -3290,15 +3316,15 @@
 
              IF ( salsa  .AND.  nest_salsa )  THEN
-                DO  ib = 1, nbins_aerosol
-                   CALL pmci_interp_1sto_lr( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),    &
+                DO  lb = 1, nbins_aerosol
+                   CALL pmci_interp_1sto_lr( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
                                              kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
                 ENDDO
-                DO  ic = 1, nbins_aerosol * ncomponents_mass
-                   CALL pmci_interp_1sto_lr( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),        &
+                DO  lc = 1, nbins_aerosol * ncomponents_mass
+                   CALL pmci_interp_1sto_lr( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
                                              kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      CALL pmci_interp_1sto_lr( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),         &
+                   DO  lg = 1, ngases_salsa
+                      CALL pmci_interp_1sto_lr( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
                                                 kcto, jflo, jfuo, kflo, kfuo, 'l', 's')
                    ENDDO
@@ -3315,11 +3341,11 @@
              CALL pmci_interp_1sto_lr( w, wc, kctw, jflo, jfuo, kflw, kfuw, 'r', 'w' )
 
-             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                        &
-                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                         &
-                     .NOT. constant_diffusion ) )  THEN
+             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
+                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
+                     .NOT.  constant_diffusion ) )  THEN
 !                CALL pmci_interp_1sto_lr( e, ec, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, 'r', 'e' )
 !
 !--             Interpolation of e is replaced by the Neumann condition. 
-                DO ibgp = nx+1, nx+nbgp
+                DO  ibgp = nx+1, nx+nbgp
                    e(nzb:nzt,nys:nyn,ibgp) = e(nzb:nzt,nys:nyn,nx)
                 ENDDO
@@ -3328,8 +3354,7 @@
              IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
                 CALL pmci_interp_1sto_lr( diss, dissc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
-
              ENDIF
 
-             IF ( .NOT. neutral )  THEN
+             IF (  .NOT.  neutral )  THEN
                 CALL pmci_interp_1sto_lr( pt, ptc, kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
              ENDIF
@@ -3362,15 +3387,15 @@
 
              IF ( salsa  .AND.  nest_salsa )  THEN
-                DO  ib = 1, nbins_aerosol
-                   CALL pmci_interp_1sto_lr( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),    &
+                DO  lb = 1, nbins_aerosol
+                   CALL pmci_interp_1sto_lr( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
                                              kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
                 ENDDO
-                DO  ic = 1, nbins_aerosol * ncomponents_mass
-                   CALL pmci_interp_1sto_lr( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),        &
+                DO  lc = 1, nbins_aerosol * ncomponents_mass
+                   CALL pmci_interp_1sto_lr( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
                                              kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      CALL pmci_interp_1sto_lr( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),         &
+                   DO  lg = 1, ngases_salsa
+                      CALL pmci_interp_1sto_lr( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
                                                 kcto, jflo, jfuo, kflo, kfuo, 'r', 's' )
                    ENDDO
@@ -3387,11 +3412,11 @@
              CALL pmci_interp_1sto_sn( u, uc, kcto, iflu, ifuu, kflo, kfuo, 's', 'u' )
 
-             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
-                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
-                     .NOT. constant_diffusion ) )  THEN
+             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       &
+                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
+                     .NOT.  constant_diffusion ) )  THEN
 !                CALL pmci_interp_1sto_sn( e, ec, kcto, iflo, ifuo, kflo, kfuo, 's', 'e' )
 !
 !--             Interpolation of e is replaced by the Neumann condition. 
-                DO jbgp = -nbgp, -1
+                DO  jbgp = -nbgp, -1
                    e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,0,nxl:nxr)
                 ENDDO
@@ -3402,5 +3427,5 @@
              ENDIF
 
-             IF ( .NOT. neutral )  THEN
+             IF (  .NOT.  neutral )  THEN
                 CALL pmci_interp_1sto_sn( pt, ptc, kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
              ENDIF
@@ -3433,15 +3458,15 @@
              
              IF ( salsa  .AND.  nest_salsa )  THEN
-                DO  ib = 1, nbins_aerosol
-                   CALL pmci_interp_1sto_sn( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),    &
+                DO  lb = 1, nbins_aerosol
+                   CALL pmci_interp_1sto_sn( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
                                              kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
                 ENDDO
-                DO  ic = 1, nbins_aerosol * ncomponents_mass
-                   CALL pmci_interp_1sto_sn( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),        &
+                DO  lc = 1, nbins_aerosol * ncomponents_mass
+                   CALL pmci_interp_1sto_sn( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
                                              kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      CALL pmci_interp_1sto_sn( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),         &
+                   DO  lg = 1, ngases_salsa
+                      CALL pmci_interp_1sto_sn( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
                                                 kcto, iflo, ifuo, kflo, kfuo, 's', 's' )
                    ENDDO
@@ -3458,11 +3483,11 @@
              CALL pmci_interp_1sto_sn( u, uc, kcto, iflu, ifuu, kflo, kfuo, 'n', 'u' )
 
-             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.                        & 
-                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                         &
-                     .NOT. constant_diffusion ) )  THEN
+             IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                       & 
+                  (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                        &
+                     .NOT.  constant_diffusion ) )  THEN
 !                CALL pmci_interp_1sto_sn( e, ec, kcto, iflo, ifuo, kflo, kfuo, 'n', 'e' )
 !
 !--             Interpolation of e is replaced by the Neumann condition. 
-                DO jbgp = ny+1, ny+nbgp
+                DO  jbgp = ny+1, ny+nbgp
                    e(nzb:nzt,jbgp,nxl:nxr) = e(nzb:nzt,ny,nxl:nxr)
                 ENDDO
@@ -3473,5 +3498,5 @@
              ENDIF
 
-             IF ( .NOT. neutral )  THEN
+             IF (  .NOT.  neutral )  THEN
                 CALL pmci_interp_1sto_sn( pt, ptc, kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
              ENDIF
@@ -3504,15 +3529,15 @@
              
              IF ( salsa  .AND.  nest_salsa )  THEN
-                DO  ib = 1, nbins_aerosol
-                   CALL pmci_interp_1sto_sn( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),    &
+                DO  lb = 1, nbins_aerosol
+                   CALL pmci_interp_1sto_sn( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),   &
                                              kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
                 ENDDO
-                DO  ic = 1, nbins_aerosol * ncomponents_mass
-                   CALL pmci_interp_1sto_sn( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),        &
+                DO  lc = 1, nbins_aerosol * ncomponents_mass
+                   CALL pmci_interp_1sto_sn( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),       &
                                              kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
                 ENDDO
                 IF ( .NOT. salsa_gases_from_chem )  THEN
-                   DO  ig = 1, ngases_salsa
-                      CALL pmci_interp_1sto_sn( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),         &
+                   DO  lg = 1, ngases_salsa
+                      CALL pmci_interp_1sto_sn( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),          &
                                                 kcto, iflo, ifuo, kflo, kfuo, 'n', 's' )
                    ENDDO
@@ -3529,7 +3554,7 @@
        CALL pmci_interp_1sto_t( v, vc, kcto, iflo, ifuo, jflv, jfuv, 'v' )
 
-       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
-            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
-               .NOT. constant_diffusion ) )  THEN
+       IF ( (         rans_mode_parent  .AND.         rans_mode )  .OR.                             &
+            (  .NOT.  rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.                              &
+               .NOT.  constant_diffusion ) )  THEN
 !          CALL pmci_interp_1sto_t( e, ec, kcto, iflo, ifuo, jflo, jfuo, 'e' )
 !
@@ -3570,15 +3595,15 @@
        
        IF ( salsa  .AND.  nest_salsa )  THEN
-          DO  ib = 1, nbins_aerosol
-             CALL pmci_interp_1sto_t( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),           &
+          DO  lb = 1, nbins_aerosol
+             CALL pmci_interp_1sto_t( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),          &
                                       kcto, iflo, ifuo, jflo, jfuo, 's' )
           ENDDO
-          DO  ic = 1, nbins_aerosol * ncomponents_mass
-             CALL pmci_interp_1sto_t( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),               &
+          DO  lc = 1, nbins_aerosol * ncomponents_mass
+             CALL pmci_interp_1sto_t( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),              &
                                       kcto, iflo, ifuo, jflo, jfuo, 's' )
           ENDDO
           IF ( .NOT. salsa_gases_from_chem )  THEN
-             DO  ig = 1, ngases_salsa
-                CALL pmci_interp_1sto_t( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),                &
+             DO  lg = 1, ngases_salsa
+                CALL pmci_interp_1sto_t( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),                 &
                                          kcto, iflo, ifuo, jflo, jfuo, 's' )
              ENDDO
@@ -3596,10 +3621,10 @@
 !--   Note that TKE is not anterpolated.
       IMPLICIT NONE
-      INTEGER(iwp) ::  ib         !< running index for aerosol size bins
-      INTEGER(iwp) ::  ic         !< running index for aerosol mass bins
-      INTEGER(iwp) ::  n          !< running index for number of chemical species
-      INTEGER(iwp) ::  ig         !< running index for salsa gases
-
-
+      INTEGER(iwp) ::  lb         !< Running index for aerosol size bins
+      INTEGER(iwp) ::  lc         !< Running index for aerosol mass bins
+      INTEGER(iwp) ::  lg         !< Running index for salsa gases
+      INTEGER(iwp) ::  n          !< Running index for number of chemical species
+
+      
       CALL pmci_anterp_tophat( u,  uc,  kcto, iflu, ifuu, jflo, jfuo, kflo, kfuo, ijkfc_u, 'u' )
       CALL pmci_anterp_tophat( v,  vc,  kcto, iflo, ifuo, jflv, jfuv, kflo, kfuo, ijkfc_v, 'v' )
@@ -3613,6 +3638,6 @@
 !--      Anterpolation of dissipation rate only if TKE-e closure is applied.
          IF ( rans_tke_e )  THEN
-            CALL pmci_anterp_tophat( diss, dissc, kcto, iflo, ifuo, jflo, jfuo,&
-                                     kflo, kfuo, ijkfc_s, 'diss' )
+            CALL pmci_anterp_tophat( diss, dissc, kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo,         &
+                 ijkfc_s, 'diss' )
          ENDIF
 
@@ -3661,15 +3686,15 @@
       
       IF ( salsa  .AND.  nest_salsa )  THEN
-         DO  ib = 1, nbins_aerosol
-            CALL pmci_anterp_tophat( aerosol_number(ib)%conc, aerosol_number_c(:,:,:,ib),            &
+         DO  lb = 1, nbins_aerosol
+            CALL pmci_anterp_tophat( aerosol_number(lb)%conc, aerosol_number_c(:,:,:,lb),           &
                                      kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
          ENDDO
-         DO  ic = 1, nbins_aerosol * ncomponents_mass
-            CALL pmci_anterp_tophat( aerosol_mass(ic)%conc, aerosol_mass_c(:,:,:,ic),                &
+         DO  lc = 1, nbins_aerosol * ncomponents_mass
+            CALL pmci_anterp_tophat( aerosol_mass(lc)%conc, aerosol_mass_c(:,:,:,lc),               &
                                      kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
          ENDDO
          IF ( .NOT. salsa_gases_from_chem )  THEN
-            DO  ig = 1, ngases_salsa
-               CALL pmci_anterp_tophat( salsa_gas(ig)%conc, salsa_gas_c(:,:,:,ig),                 &
+            DO  lg = 1, ngases_salsa
+               CALL pmci_anterp_tophat( salsa_gas(lg)%conc, salsa_gas_c(:,:,:,lg),                  &
                                         kcto, iflo, ifuo, jflo, jfuo, kflo, kfuo, ijkfc_s, 's' )
             ENDDO
@@ -3681,5 +3706,5 @@
 
 
-   SUBROUTINE pmci_interp_1sto_lr( f, fc, kct, jfl, jfu, kfl, kfu, edge, var )
+   SUBROUTINE pmci_interp_1sto_lr( child_array, parent_array, kct, jfl, jfu, kfl, kfu, edge, var )
 !
 !--   Interpolation of ghost-node values used as the child-domain boundary
@@ -3687,30 +3712,42 @@
       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
-      INTEGER(iwp) :: kct                                     !< The parent-grid index in z-direction just below the boundary value node
-      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 end 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 end index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
+      
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - y direction
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - y direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - z direction
+
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
+      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
+
       CHARACTER(LEN=1), INTENT(IN) ::  edge                   !< Edge symbol: 'l' or 'r'
-      CHARACTER(LEN=1), INTENT(IN) ::  var                    !< Variable symbol: 'u', 'v', 'w' or 's'
+      CHARACTER(LEN=1), INTENT(IN) ::  var                    !< Variable symbol: 'u', 'v', 'w' or 's'      
 !
 !--   Local variables:
-      INTEGER(iwp) ::  ib       !< Fixed i-index pointing to the node just behind the boundary-value node
-      INTEGER(iwp) ::  ibc      !< Fixed i-index pointing to the boundary-value nodes (either i or iend)
-      INTEGER(iwp) ::  ibgp     !< Index running over the redundant boundary ghost points in i-direction
+      INTEGER(iwp) ::  icb      !< Fixed child-grid index in the x-direction pointing to the node just behind the
+                                !< boundary-value node
+      INTEGER(iwp) ::  icbc     !< Fixed child-grid index in the x-direction pointing to the boundary-value nodes
+      INTEGER(iwp) ::  icbgp    !< Index running over the redundant boundary ghost points in the x-direction
       INTEGER(iwp) ::  ierr     !< MPI error code
-      INTEGER(iwp) ::  j        !< Running index in the y-direction
-      INTEGER(iwp) ::  k        !< Running index in the z-direction
-      INTEGER(iwp) ::  lbeg     !< l-index pointing to the starting point of workarrc_lr in the l-direction 
-      INTEGER(iwp) ::  lw       !< Reduced l-index for workarrc_lr pointing to the boundary ghost node
-      INTEGER(iwp) ::  lwp      !< Reduced l-index for workarrc_lr pointing to the first prognostic node
-      INTEGER(iwp) ::  m        !< Parent-grid running index in the y-direction
-      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  ipbeg    !< Parent-grid index in the x-direction pointing to the starting point of workarr_lr
+                                !< in the parent-grid array
+      INTEGER(iwp) ::  ipw      !< Reduced parent-grid index in the x-direction for workarr_lr pointing to
+                                !< the boundary ghost node
+      INTEGER(iwp) ::  ipwp     !< Reduced parent-grid index in the x-direction for workarr_lr pointing to
+                                !< the first prognostic node
+      INTEGER(iwp) ::  jc       !< Running child-grid index in the y-direction
+      INTEGER(iwp) ::  jp       !< Running parent-grid index in the y-direction
+      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction
+      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction      
+      
       REAL(wp) ::  cb           !< Interpolation coefficient for the boundary ghost node  
       REAL(wp) ::  cp           !< Interpolation coefficient for the first prognostic node
-      REAL(wp) ::  f_interp_1   !< Value interpolated in x direction from the parent-grid data
-      REAL(wp) ::  f_interp_2   !< Auxiliary value interpolated in x direction from the parent-grid data
+      REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
+      REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
 
 ! 
@@ -3720,55 +3757,53 @@
 !--      For u, nxl is a ghost node, but not for the other variables
          IF ( var == 'u' )  THEN
-            ibc   = nxl    
-            ib    = ibc - 1
-            lw    = 2
-            lwp   = lw        ! This is redundant when var == 'u' 
-            lbeg  = icl
+            icbc  = nxl    
+            icb   = icbc - 1
+            ipw   = 2
+            ipwp  = ipw        ! This is redundant when var == 'u' 
+            ipbeg = ipl
          ELSE
-            ibc   = nxl - 1
-            ib    = ibc - 1
-            lw    = 1
-            lwp   = 2
-            lbeg  = icl
+            icbc  = nxl - 1
+            icb   = icbc - 1
+            ipw   = 1
+            ipwp  = 2
+            ipbeg = ipl
          ENDIF        
       ELSEIF ( edge == 'r' )  THEN
          IF ( var == 'u' )  THEN
-            ibc   = nxr + 1 
-            ib    = ibc + 1
-            lw    = 1
-            lwp   = lw        ! This is redundant when var == 'u'            
-            lbeg  = icr - 2
+            icbc  = nxr + 1 
+            icb   = icbc + 1
+            ipw   = 1
+            ipwp  = ipw        ! This is redundant when var == 'u'            
+            ipbeg = ipr - 2
          ELSE
-            ibc   = nxr + 1
-            ib    = ibc + 1
-            lw    = 1
-            lwp   = 0
-            lbeg  = icr - 2
+            icbc  = nxr + 1
+            icb   = icbc + 1
+            ipw   = 1
+            ipwp  = 0
+            ipbeg = ipr - 2
          ENDIF         
       ENDIF
 !
 !--   Interpolation coefficients
-      IF  ( interpolation_scheme_lrsn == 1 )  THEN
+      IF ( interpolation_scheme_lrsn == 1 )  THEN
          cb = 1.0_wp  ! 1st-order upwind 
-      ELSE IF  ( interpolation_scheme_lrsn == 2 )  THEN
+      ELSE IF ( interpolation_scheme_lrsn == 2 )  THEN
          cb = 0.5_wp  ! 2nd-order central
       ELSE
          cb = 0.5_wp  ! 2nd-order central (default)
       ENDIF         
-      cp    = 1.0_wp - cb
-!
-!--   Substitute the necessary parent-grid data to the work array workarrc_lr. 
-      workarrc_lr = 0.0_wp     
-      IF  ( pdims(2) > 1 )  THEN
+      cp = 1.0_wp - cb
+!
+!--   Substitute the necessary parent-grid data to the work array workarr_lr. 
+      workarr_lr = 0.0_wp     
+      IF ( pdims(2) > 1 )  THEN
 #if defined( __parallel )
-         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)
+         IF ( bc_dirichlet_s )  THEN 
+            workarr_lr(0:pg%nz+1,jpsw:jpnw-1,0:2) = parent_array(0:pg%nz+1,jpsw:jpnw-1,ipbeg:ipbeg+2)
+         ELSE IF ( bc_dirichlet_n )  THEN
+            workarr_lr(0:pg%nz+1,jpsw+1:jpnw,0:2) = parent_array(0:pg%nz+1,jpsw+1:jpnw,ipbeg:ipbeg+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)
+            workarr_lr(0:pg%nz+1,jpsw+1:jpnw-1,0:2)                                                 &
+                 = parent_array(0:pg%nz+1,jpsw+1:jpnw-1,ipbeg:ipbeg+2)
          ENDIF
 !
@@ -3778,30 +3813,25 @@
 !--      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)            
+         CALL MPI_SENDRECV( workarr_lr(0,jpsw+1,0), 1, workarr_lr_exchange_type, psouth,  0,        &
+                            workarr_lr(0,jpnw,0), 1, workarr_lr_exchange_type, pnorth,  0, comm2d,  &
+                            status, ierr )                              
+!                                                                             
+!--      From north to south                                                  
+         CALL MPI_SENDRECV( workarr_lr(0,jpnw-1,0), 1, workarr_lr_exchange_type, pnorth,  1,        &
+                            workarr_lr(0,jpsw,0), 1, workarr_lr_exchange_type, psouth,  1, comm2d,  &
+                            status, ierr )                               
+#endif                                                                        
+      ELSE                                                                    
+         workarr_lr(0:pg%nz+1,jpsw:jpnw,0:2) = parent_array(0:pg%nz+1,jpsw:jpnw,ipbeg:ipbeg+2)            
       ENDIF
 
-      IF  ( var == 'u' )  THEN
-
-         DO  m = jcsw, jcnw
-            DO n = 0, kct 
+      IF ( var == 'u' )  THEN
+
+         DO  jp = jpsw, jpnw
+            DO  kp = 0, kct 
                
-               DO  j = jfl(m), jfu(m)
-                  DO  k = kfl(n), kfu(n)
-                     f(k,j,ibc) = workarrc_lr(n,m,lw)
+               DO  jc = jfl(jp), jfu(jp)
+                  DO  kc = kfl(kp), kfu(kp)
+                     child_array(kc,jc,icbc) = workarr_lr(kp,jp,ipw)
                   ENDDO
                ENDDO
@@ -3810,21 +3840,21 @@
          ENDDO
 
-      ELSE IF  ( var == 'v' )  THEN
+      ELSE IF ( var == 'v' )  THEN
          
-         DO  m = jcsw, jcnw-1
-            DO n = 0, kct 
+         DO  jp = jpsw, jpnw-1
+            DO  kp = 0, kct 
 !
 !--            First interpolate to the flux point
-               f_interp_1 = cb * workarrc_lr(n,m,lw)   + cp * workarrc_lr(n,m,lwp)
-               f_interp_2 = cb * workarrc_lr(n,m+1,lw) + cp * workarrc_lr(n,m+1,lwp)
+               c_interp_1 = cb * workarr_lr(kp,jp,ipw)   + cp * workarr_lr(kp,jp,ipwp)
+               c_interp_2 = cb * workarr_lr(kp,jp+1,ipw) + cp * workarr_lr(kp,jp+1,ipwp)
 !
 !--            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 * ( f_interp_1 + f_interp_2 )
+               DO  jc = jfl(jp), jfl(jp+1)
+                  child_array(kfl(kp):kfu(kp),jc,icbc) = 0.5_wp * ( c_interp_1 + c_interp_2 )
                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) = f_interp_1
+               IF  ( MOD( jfl(jp), jgsr ) == 0 )  THEN
+                  child_array(kfl(kp):kfu(kp),jfl(jp),icbc) = c_interp_1
                ENDIF
             ENDDO
@@ -3832,8 +3862,8 @@
 !
 !--      Finally, set the values along the last matching grid-line  
-         IF  ( MOD( jfl(jcnw), jgsr ) == 0 )  THEN
-            DO  n = 0, kct
-               f_interp_1 = cb * workarrc_lr(n,jcnw,lw) + cp * workarrc_lr(n,jcnw,lwp)
-               f(kfl(n):kfu(n),jfl(jcnw),ibc) = f_interp_1
+         IF ( MOD( jfl(jpnw), jgsr ) == 0 )  THEN
+            DO  kp = 0, kct
+               c_interp_1 = cb * workarr_lr(kp,jpnw,ipw) + cp * workarr_lr(kp,jpnw,ipwp)
+               child_array(kfl(kp):kfu(kp),jfl(jpnw),icbc) = c_interp_1
             ENDDO
          ENDIF
@@ -3843,33 +3873,34 @@
 !--      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)
+         IF ( jfl(jpnw) < nyn )  THEN
+            DO  kp = 0, kct
+               DO  jc = jfl(jpnw) + 1, nyn
+                  child_array(kfl(kp):kfu(kp),jc,icbc) = child_array(kfl(kp):kfu(kp),jfl(jpnw),icbc)
                ENDDO
             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  
+      ELSE IF ( var == 'w' )  THEN
+
+         DO  jp = jpsw, jpnw
+            DO  kp = 0, kct + 1   ! It is important to go up to kct+1  
 !
 !--            Interpolate to the flux point
-               f_interp_1 = cb * workarrc_lr(n,m,lw) + cp * workarrc_lr(n,m,lwp)
+               c_interp_1 = cb * workarr_lr(kp,jp,ipw) + cp * workarr_lr(kp,jp,ipwp)
 !
 !--            First substitute only the matching-node values
-               f(kfu(n),jfl(m):jfu(m),ibc) = f_interp_1
+               child_array(kfu(kp),jfl(jp):jfu(jp),icbc) = c_interp_1
                   
             ENDDO
          ENDDO
 
-         DO  m = jcsw, jcnw
-            DO n = 1, kct + 1   ! It is important to go up to kct+1  
+         DO  jp = jpsw, jpnw
+            DO  kp = 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) )
+               DO  kc = kfu(kp-1) + 1, kfu(kp) - 1 
+                  child_array(kc,jfl(jp):jfu(jp),icbc) =                                            &
+                       0.5_wp * ( child_array(kfu(kp-1),jfl(jp):jfu(jp),icbc)                       &
+                       + child_array(kfu(kp),jfl(jp):jfu(jp),icbc) )
                ENDDO
                   
@@ -3879,12 +3910,12 @@
       ELSE   ! any scalar
          
-         DO  m = jcsw, jcnw
-            DO n = 0, kct 
+         DO  jp = jpsw, jpnw
+            DO  kp = 0, kct 
 !
 !--            Interpolate to the flux point
-               f_interp_1 = cb * workarrc_lr(n,m,lw) + cp * workarrc_lr(n,m,lwp)
-               DO  j = jfl(m), jfu(m)
-                  DO  k = kfl(n), kfu(n)
-                     f(k,j,ibc) = f_interp_1
+               c_interp_1 = cb * workarr_lr(kp,jp,ipw) + cp * workarr_lr(kp,jp,ipwp)
+               DO  jc = jfl(jp), jfu(jp)
+                  DO  kc = kfl(kp), kfu(kp)
+                     child_array(kc,jc,icbc) = c_interp_1
                   ENDDO
                ENDDO
@@ -3897,10 +3928,10 @@
 !--   Fill up also the redundant 2nd and 3rd ghost-node layers
       IF ( edge == 'l' )  THEN
-         DO  ibgp = -nbgp, ib
-            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,ibc)
+         DO  icbgp = -nbgp, icb
+            child_array(0:nzt+1,nysg:nyng,icbgp) = child_array(0:nzt+1,nysg:nyng,icbc)
          ENDDO
       ELSEIF ( edge == 'r' )  THEN
-         DO  ibgp = ib, nx+nbgp
-            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,ibc)
+         DO  icbgp = icb, nx+nbgp
+            child_array(0:nzt+1,nysg:nyng,icbgp) = child_array(0:nzt+1,nysg:nyng,icbc)
          ENDDO
       ENDIF
@@ -3910,5 +3941,5 @@
 
 
-   SUBROUTINE pmci_interp_1sto_sn( f, fc, kct, ifl, ifu, kfl, kfu, edge, var )
+   SUBROUTINE pmci_interp_1sto_sn( child_array, parent_array, kct, ifl, ifu, kfl, kfu, edge, var )
 !
 !--   Interpolation of ghost-node values used as the child-domain boundary
@@ -3916,31 +3947,42 @@
       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
-      INTEGER(iwp) :: kct                                     !< The parent-grid index in z-direction just below the boundary value node
-      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(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 end index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
+      
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - z direction
+      
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
+      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
+
       CHARACTER(LEN=1), INTENT(IN) ::  edge   !< Edge symbol: 's' or 'n' 
       CHARACTER(LEN=1), INTENT(IN) ::  var    !< Variable symbol: 'u', 'v', 'w' or 's'
 !
 !--   Local variables:      
-      INTEGER(iwp) ::  i        !< Running index in the x-direction
+      INTEGER(iwp) ::  ic       !< Running child-grid index in the x-direction
       INTEGER(iwp) ::  ierr     !< MPI error code
-      INTEGER(iwp) ::  jb       !< Fixed j-index pointing to the node just behind the boundary-value node
-      INTEGER(iwp) ::  jbc      !< Fixed j-index pointing to the boundary-value nodes (either j or jend)
-      INTEGER(iwp) ::  jbgp     !< Index running over the redundant boundary ghost points in j-direction
-      INTEGER(iwp) ::  k        !< Running index in the z-direction
-      INTEGER(iwp) ::  l        !< Parent-grid running index in the x-direction
-      INTEGER(iwp) ::  mbeg     !< m-index pointing to the starting point of workarrc_sn in the m-direction 
-      INTEGER(iwp) ::  mw       !< Reduced m-index for workarrc_sn pointing to the boundary ghost node 
-      INTEGER(iwp) ::  mwp      !< Reduced m-index for workarrc_sn pointing to the first prognostic node
-      INTEGER(iwp) ::  n        !< Parent-grid running index in the z-direction
+      INTEGER(iwp) ::  ip       !< Running parent-grid index in the x-direction
+      INTEGER(iwp) ::  jcb      !< Fixed child-grid index in the y-direction pointing to the node just behind the
+                                !< boundary-value node
+      INTEGER(iwp) ::  jcbc     !< Fixed child-grid index in the y-direction pointing to the boundary-value nodes
+      INTEGER(iwp) ::  jcbgp    !< Index running over the redundant boundary ghost points in y-direction
+      INTEGER(iwp) ::  kc       !< Running child-grid index in the z-direction      
+      INTEGER(iwp) ::  kp       !< Running parent-grid index in the z-direction
+      INTEGER(iwp) ::  jpbeg    !< Parent-grid index in the y-direction pointing to the starting point of workarr_sn
+                                !< in the parent-grid array
+      INTEGER(iwp) ::  jpw      !< Reduced parent-grid index in the y-direction for workarr_sn pointing to
+                                !< the boundary ghost node 
+      INTEGER(iwp) ::  jpwp     !< Reduced parent-grid index in the y-direction for workarr_sn pointing to
+                                !< the first prognostic node
       REAL(wp) ::  cb           !< Interpolation coefficient for the boundary ghost node  
       REAL(wp) ::  cp           !< Interpolation coefficient for the first prognostic node
-      REAL(wp) ::  f_interp_1   !< Value interpolated in y direction from the parent-grid data
-      REAL(wp) ::  f_interp_2   !< Auxiliary value interpolated in y direction from the parent-grid data
-
+      REAL(wp) ::  c_interp_1   !< Value interpolated to the flux point in x direction from the parent-grid data
+      REAL(wp) ::  c_interp_2   !< Auxiliary value interpolated  to the flux point in x direction from the parent-grid data
+      
 !
 !--   Check which edge is to be handled: south or north
@@ -3949,55 +3991,53 @@
 !--      For v, nys is a ghost node, but not for the other variables
          IF ( var == 'v' )  THEN
-            jbc   = nys
-            jb    = jbc - 1
-            mw    = 2
-            mwp   = 2         ! This is redundant when var == 'v'
-            mbeg  = jcs
+            jcbc  = nys
+            jcb   = jcbc - 1
+            jpw   = 2
+            jpwp  = 2         ! This is redundant when var == 'v'
+            jpbeg = jps
          ELSE
-            jbc   = nys - 1
-            jb    = jbc - 1
-            mw    = 1
-            mwp   = 2
-            mbeg  = jcs
+            jcbc  = nys - 1
+            jcb   = jcbc - 1
+            jpw   = 1
+            jpwp  = 2
+            jpbeg = jps
          ENDIF
       ELSEIF ( edge == 'n' )  THEN
          IF ( var == 'v' )  THEN
-            jbc   = nyn + 1
-            jb    = jbc + 1
-            mw    = 1
-            mwp   = 0         ! This is redundant when var == 'v'     
-            mbeg  = jcn - 2
+            jcbc  = nyn + 1
+            jcb   = jcbc + 1
+            jpw   = 1
+            jpwp  = 0         ! This is redundant when var == 'v'     
+            jpbeg = jpn - 2
          ELSE
-            jbc   = nyn + 1
-            jb    = jbc + 1
-            mw    = 1
-            mwp   = 0
-            mbeg  = jcn - 2
+            jcbc  = nyn + 1
+            jcb   = jcbc + 1
+            jpw   = 1
+            jpwp  = 0
+            jpbeg = jpn - 2
          ENDIF
       ENDIF
 !
 !--   Interpolation coefficients
-      IF  ( interpolation_scheme_lrsn == 1 )  THEN
+      IF ( interpolation_scheme_lrsn == 1 )  THEN
          cb = 1.0_wp  ! 1st-order upwind 
-      ELSE IF  ( interpolation_scheme_lrsn == 2 )  THEN
+      ELSE IF ( interpolation_scheme_lrsn == 2 )  THEN
          cb = 0.5_wp  ! 2nd-order central
       ELSE
          cb = 0.5_wp  ! 2nd-order central (default)
       ENDIF         
-      cp    = 1.0_wp - cb
-!
-!--   Substitute the necessary parent-grid data to the work array workarrc_sn. 
-      workarrc_sn = 0.0_wp     
-      IF  ( pdims(1) > 1 )  THEN
+      cp = 1.0_wp - cb
+!
+!--   Substitute the necessary parent-grid data to the work array workarr_sn. 
+      workarr_sn = 0.0_wp     
+      IF ( pdims(1) > 1 )  THEN
 #if defined( __parallel )
-         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)
+         IF ( bc_dirichlet_l )  THEN
+            workarr_sn(0:pg%nz+1,0:2,iplw:iprw-1) = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw:iprw-1)
+         ELSE IF ( bc_dirichlet_r )  THEN
+            workarr_sn(0:pg%nz+1,0:2,iplw+1:iprw) = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw+1:iprw)
          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)
+            workarr_sn(0:pg%nz+1,0:2,iplw+1:iprw-1)                                                 &
+                 = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw+1:iprw-1)
          ENDIF
 !
@@ -4007,30 +4047,26 @@
 !--      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)
+         CALL MPI_SENDRECV( workarr_sn(0,0,iplw+1), 1, workarr_sn_exchange_type, pleft,   0,        &
+                            workarr_sn(0,0,iprw), 1, workarr_sn_exchange_type, pright, 0, comm2d,   &
+                            status, ierr )
+!                                                                           
+!--      From right to left                                                 
+         CALL MPI_SENDRECV( workarr_sn(0,0,iprw-1), 1, workarr_sn_exchange_type, pright,  1,        &
+                            workarr_sn(0,0,iplw), 1, workarr_sn_exchange_type, pleft, 1, comm2d,    &
+                            status, ierr )
+#endif                                                                      
+      ELSE                                                                  
+         workarr_sn(0:pg%nz+1,0:2,iplw+1:iprw-1)                                                    &
+              = parent_array(0:pg%nz+1,jpbeg:jpbeg+2,iplw+1:iprw-1)
       ENDIF
 
-      IF  ( var == 'v' )  THEN
-
-         DO  l = iclw, icrw
-            DO n = 0, kct 
+      IF ( var == 'v' )  THEN
+
+         DO  ip = iplw, iprw
+            DO  kp = 0, kct 
                
-               DO  i = ifl(l), ifu(l)
-                  DO  k = kfl(n), kfu(n)
-                     f(k,jbc,i) = workarrc_sn(n,mw,l)
+               DO  ic = ifl(ip), ifu(ip)
+                  DO  kc = kfl(kp), kfu(kp)
+                     child_array(kc,jcbc,ic) = workarr_sn(kp,jpw,ip)
                   ENDDO
                ENDDO
@@ -4039,21 +4075,21 @@
          ENDDO
 
-      ELSE IF  ( var == 'u' )  THEN
+      ELSE IF ( var == 'u' )  THEN
          
-         DO  l = iclw, icrw-1
-            DO n = 0, kct
+         DO  ip = iplw, iprw - 1
+            DO  kp = 0, kct
 !
 !--            First interpolate to the flux point
-               f_interp_1 = cb * workarrc_sn(n,mw,l)   + cp * workarrc_sn(n,mwp,l)
-               f_interp_2 = cb * workarrc_sn(n,mw,l+1) + cp * workarrc_sn(n,mwp,l+1)
+               c_interp_1 = cb * workarr_sn(kp,jpw,ip)   + cp * workarr_sn(kp,jpwp,ip)
+               c_interp_2 = cb * workarr_sn(kp,jpw,ip+1) + cp * workarr_sn(kp,jpwp,ip+1)
 !
 !--            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 * ( f_interp_1 + f_interp_2 )
+               DO  ic = ifl(ip), ifl(ip+1)
+                  child_array(kfl(kp):kfu(kp),jcbc,ic) = 0.5_wp * ( c_interp_1 + c_interp_2 )
                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)) = f_interp_1
+               IF ( MOD( ifl(ip), igsr ) == 0 )  THEN
+                  child_array(kfl(kp):kfu(kp),jcbc,ifl(ip)) = c_interp_1
                ENDIF
 
@@ -4062,8 +4098,8 @@
 !
 !--      Finally, set the values along the last matching grid-line  
-         IF  ( MOD( ifl(icrw), igsr ) == 0 )  THEN
-            DO  n = 0, kct
-               f_interp_1 = cb * workarrc_sn(n,mw,icrw) + cp * workarrc_sn(n,mwp,icrw)
-               f(kfl(n):kfu(n),jbc,ifl(icrw)) = f_interp_1
+         IF ( MOD( ifl(iprw), igsr ) == 0 )  THEN
+            DO  kp = 0, kct
+               c_interp_1 = cb * workarr_sn(kp,jpw,iprw) + cp * workarr_sn(kp,jpwp,iprw)
+               child_array(kfl(kp):kfu(kp),jcbc,ifl(iprw)) = c_interp_1
             ENDDO
          ENDIF
@@ -4073,33 +4109,34 @@
 !--      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))
+         IF ( ifl(iprw) < nxr )  THEN
+            DO  kp = 0, kct
+               DO  ic = ifl(iprw) + 1, nxr
+                  child_array(kfl(kp):kfu(kp),jcbc,ic) = child_array(kfl(kp):kfu(kp),jcbc,ifl(iprw))
                ENDDO
             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  
+      ELSE IF ( var == 'w' )  THEN
+
+         DO  ip = iplw, iprw
+            DO  kp = 0, kct + 1   ! It is important to go up to kct+1  
 !
 !--            Interpolate to the flux point
-               f_interp_1 = cb * workarrc_sn(n,mw,l) + cp * workarrc_sn(n,mwp,l)
-!                  
-!--            First substitute only the matching-node values                 
-               f(kfu(n),jbc,ifl(l):ifu(l)) = f_interp_1
+               c_interp_1 = cb * workarr_sn(kp,jpw,ip) + cp * workarr_sn(kp,jpwp,ip)
+!
+!--            First substitute only the matching-node values
+               child_array(kfu(kp),jcbc,ifl(ip):ifu(ip)) = c_interp_1
 
             ENDDO
          ENDDO
 
-         DO  l = iclw, icrw
-            DO n = 1, kct + 1   ! It is important to go up to kct+1  
+         DO  ip = iplw, iprw
+            DO  kp = 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)) )
+               DO  kc = kfu(kp-1) + 1, kfu(kp) - 1 
+                  child_array(kc,jcbc,ifl(ip):ifu(ip)) =                                            &
+                       0.5_wp * ( child_array(kfu(kp-1),jcbc,ifl(ip):ifu(ip))                       &
+                       + child_array(kfu(kp),jcbc,ifl(ip):ifu(ip)) )
                ENDDO
 
@@ -4109,12 +4146,12 @@
       ELSE   ! Any scalar
          
-         DO  l = iclw, icrw
-            DO n = 0, kct 
+         DO  ip = iplw, iprw
+            DO  kp = 0, kct 
 !
 !--            Interpolate to the flux point
-               f_interp_1 = cb * workarrc_sn(n,mw,l) + cp * workarrc_sn(n,mwp,l)
-               DO  i = ifl(l), ifu(l)
-                  DO  k = kfl(n), kfu(n)
-                     f(k,jbc,i) = f_interp_1
+               c_interp_1 = cb * workarr_sn(kp,jpw,ip) + cp * workarr_sn(kp,jpwp,ip)
+               DO  ic = ifl(ip), ifu(ip)
+                  DO  kc = kfl(kp), kfu(kp)
+                     child_array(kc,jcbc,ic) = c_interp_1
                   ENDDO
                ENDDO
@@ -4127,10 +4164,10 @@
 !--   Fill up also the redundant 2nd and 3rd ghost-node layers
       IF ( edge == 's' )  THEN
-         DO  jbgp = -nbgp, jb
-            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,jbc,nxlg:nxrg)
+         DO  jcbgp = -nbgp, jcb
+            child_array(0:nzt+1,jcbgp,nxlg:nxrg) = child_array(0:nzt+1,jcbc,nxlg:nxrg)
          ENDDO
       ELSEIF ( edge == 'n' )  THEN
-         DO  jbgp = jb, ny+nbgp
-            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,jbc,nxlg:nxrg)
+         DO  jcbgp = jcb, ny+nbgp
+            child_array(0:nzt+1,jcbgp,nxlg:nxrg) = child_array(0:nzt+1,jcbc,nxlg:nxrg)
          ENDDO
       ENDIF
@@ -4140,5 +4177,5 @@
 
 
-   SUBROUTINE pmci_interp_1sto_t( f, fc, kct, ifl, ifu, jfl, jfu, var )
+   SUBROUTINE pmci_interp_1sto_t( child_array, parent_array, kct, ifl, ifu, jfl, jfu, var )
 !
 !--   Interpolation of ghost-node values used as the child-domain boundary
@@ -4146,45 +4183,57 @@
       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
-      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 end index of child cells belonging to certain parent cell - y direction
-      INTEGER(iwp) :: kct                                        !< The parent-grid index in z-direction just below the boundary value node
-      CHARACTER(LEN=1), INTENT(IN) :: var                        !< Variable symbol: 'u', 'v', 'w' or 's'
-!
-!--   Local variables:
+      INTEGER(iwp), INTENT(IN) ::  kct  !< The parent-grid index in z-direction just below the boundary value node
+      
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - y direction
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - y direction
+
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  child_array   !< Child-grid array
+      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(IN)        ::  parent_array  !< Parent-grid array
+
+      CHARACTER(LEN=1), INTENT(IN) ::  var                    !< Variable symbol: 'u', 'v', 'w' or 's'
+!
+!--   Local variables:      
+      INTEGER(iwp) ::  ic          !< Running child-grid index in the x-direction
+      INTEGER(iwp) ::  ierr        !< MPI error code
+      INTEGER(iwp) ::  iplc        !< Lower parent-grid index limit in the x-direction for copying parent-grid
+                                   !< array data to workarr_t
+      INTEGER(iwp) ::  iprc        !< Upper parent-grid index limit in the x-direction for copying parent-grid
+                                   !< array data to workarr_t
+      INTEGER(iwp) ::  jc          !< Running child-grid index in the y-direction
+      INTEGER(iwp) ::  jpsc        !< Lower parent-grid index limit in the y-direction for copying parent-grid
+                                   !< array data to workarr_t
+      INTEGER(iwp) ::  jpnc        !< Upper parent-grid-index limit in the y-direction for copying parent-grid
+                                   !< array data to workarr_t
+      INTEGER(iwp) ::  kc          !< Vertical child-grid index fixed to the boundary-value level 
+      INTEGER(iwp) ::  ip          !< Running parent-grid index in the x-direction
+      INTEGER(iwp) ::  jp          !< Running parent-grid index in the y-direction
+      INTEGER(iwp) ::  kpw         !< Reduced parent-grid index in the z-direction for workarr_t pointing to
+                                   !< the boundary ghost node
       REAL(wp)     ::  c31         !< Interpolation coefficient for the 3rd-order WS scheme
       REAL(wp)     ::  c32         !< Interpolation coefficient for the 3rd-order WS scheme
       REAL(wp)     ::  c33         !< Interpolation coefficient for the 3rd-order WS scheme
-      REAL(wp)     ::  f_interp_1  !< Value interpolated in z direction from the parent-grid data
-      REAL(wp)     ::  f_interp_2  !< Auxiliary value interpolated in z direction from the parent-grid data
-      INTEGER(iwp) ::  i           !< Child-grid index in x-direction
-      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) ::  ierr        !< MPI error code
-      INTEGER(iwp) ::  j           !< Child-grid index in y-direction
-      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) ::  k           !< Vertical child-grid index fixed to the boundary-value level 
-      INTEGER(iwp) ::  l           !< Parent-grid index in x-direction
-      INTEGER(iwp) ::  m           !< Parent-grid index in y-direction
-      INTEGER(iwp) ::  nw          !< Reduced n-index for workarrc_t pointing to the boundary ghost node 
+      REAL(wp)     ::  c_interp_1  !< Value interpolated to the flux point in z direction from the parent-grid data
+      REAL(wp)     ::  c_interp_2  !< Auxiliary value interpolated to the flux point in z direction from the parent-grid data
 
 
       IF ( var == 'w' )  THEN
-         k    = nzt
+         kc = nzt
       ELSE
-         k    = nzt + 1
+         kc = nzt + 1
       ENDIF
-      nw = 1
+      kpw = 1
 !
 !--   Interpolation coefficients
-      IF  ( interpolation_scheme_t == 1 )  THEN
+      IF ( interpolation_scheme_t == 1 )  THEN
          c31 =  0.0_wp           ! 1st-order upwind
          c32 =  1.0_wp
          c33 =  0.0_wp
-      ELSE IF  ( interpolation_scheme_t == 2 )  THEN
+      ELSE IF ( interpolation_scheme_t == 2 )  THEN
          c31 =  0.5_wp           ! 2nd-order central
          c32 =  0.5_wp
@@ -4197,26 +4246,25 @@
 !
 !--   Substitute the necessary parent-grid data to the work array. 
-!--   Note that the dimension of workarrc_t is (0:2,jcsw:jcnw,iclw:icrw),
+!--   Note that the dimension of workarr_t is (0:2,jpsw:jpnw,iplw:iprw),
 !--   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
+      iplc = iplw + 1
+      iprc = iprw - 1      
+      jpsc = jpsw + 1
+      jpnc = jpnw - 1
+      IF ( bc_dirichlet_l )  THEN
+         iplc = iplw
       ENDIF
-      IF  ( bc_dirichlet_r )  THEN
-         icrc = icrw
+      IF ( bc_dirichlet_r )  THEN
+         iprc = iprw
       ENDIF
-      IF  ( bc_dirichlet_s )  THEN
-         jcsc = jcsw
+      IF ( bc_dirichlet_s )  THEN
+         jpsc = jpsw
       ENDIF
-      IF  ( bc_dirichlet_n )  THEN
-         jcnc = jcnw
+      IF ( bc_dirichlet_n )  THEN
+         jpnc = jpnw
       ENDIF
-      workarrc_t = 0.0_wp
-!      workarrc_t(-2:3,jcsc:jcnc,iclc:icrc) = fc(kct-2:kct+3,jcsc:jcnc,iclc:icrc)
-      workarrc_t(0:2,jcsc:jcnc,iclc:icrc) = fc(kct:kct+2,jcsc:jcnc,iclc:icrc)
+      workarr_t = 0.0_wp
+      workarr_t(0:2,jpsc:jpnc,iplc:iprc) = parent_array(kct:kct+2,jpsc:jpnc,iplc:iprc)
 !
 !--   Left-right exchange if more than one PE subdomain in the x-direction.
@@ -4224,49 +4272,41 @@
 !--   not exchanged because the nest domain is not cyclic.
 #if defined( __parallel )
-      IF  ( pdims(1) > 1 )  THEN
+      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 )
+         CALL MPI_SENDRECV( workarr_t(0,jpsw,iplw+1), 1, workarr_t_exchange_type_y, pleft, 0,       &
+                            workarr_t(0,jpsw,iprw), 1, workarr_t_exchange_type_y, pright, 0,        &
+                            comm2d, status, ierr )
+!                                                                              
+!--      From right to left                                                    
+         CALL MPI_SENDRECV( workarr_t(0,jpsw,iprw-1), 1, workarr_t_exchange_type_y, pright, 1,      &
+                            workarr_t(0,jpsw,iplw), 1, workarr_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( workarr_t(0,jpsw+1,iplw), 1, workarr_t_exchange_type_x, psouth, 2,      &
+                            workarr_t(0,jpnw,iplw), 1, workarr_t_exchange_type_x, pnorth, 2,        &
+                            comm2d, status, ierr )                                           
+!                                                                              
+!--      From north to south                                                   
+         CALL MPI_SENDRECV( workarr_t(0,jpnw-1,iplw), 1, workarr_t_exchange_type_x, pnorth, 3,      &
+                            workarr_t(0,jpsw,iplw), 1, workarr_t_exchange_type_x, psouth, 3,        &
+                            comm2d, status, ierr )
       ENDIF
 #endif      
 
       IF  ( var == 'w' )  THEN
-         DO  l = iclw, icrw
-            DO  m = jcsw, jcnw
+         DO  ip = iplw, iprw
+            DO  jp = jpsw, jpnw
  
-               DO  i = ifl(l), ifu(l)
-                  DO  j = jfl(m), jfu(m)
-                     f(k,j,i) = workarrc_t(nw,m,l)
+               DO  ic = ifl(ip), ifu(ip)
+                  DO  jc = jfl(jp), jfu(jp)
+                     child_array(kc,jc,ic) = workarr_t(kpw,jp,ip)
                   ENDDO
                ENDDO
@@ -4277,25 +4317,25 @@
       ELSE IF  ( var == 'u' )  THEN
 
-         DO  l = iclw, icrw-1
-            DO  m = jcsw, jcnw
+         DO  ip = iplw, iprw - 1
+            DO  jp = jpsw, jpnw
 !
 !--            First interpolate to the flux point using the 3rd-order WS scheme
-               f_interp_1 = c31 * workarrc_t(nw-1,m,l)   + c32 * workarrc_t(nw,m,l)   + c33 * workarrc_t(nw+1,m,l)
-               f_interp_2 = c31 * workarrc_t(nw-1,m,l+1) + c32 * workarrc_t(nw,m,l+1) + c33 * workarrc_t(nw+1,m,l+1)
+               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip)   + c32 * workarr_t(kpw,jp,ip)             &
+                          + c33 * workarr_t(kpw+1,jp,ip)
+               c_interp_2 = c31 * workarr_t(kpw-1,jp,ip+1) + c32 * workarr_t(kpw,jp,ip+1)           &
+                          + c33 * workarr_t(kpw+1,jp,ip+1)
 !
 !--            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) )
-                  f(k,jfl(m):jfu(m),i) = 0.5_wp * ( f_interp_1 + f_interp_2 )
+               DO  ic = ifl(ip), ifl(ip+1)
+                  child_array(kc,jfl(jp):jfu(jp),ic) = 0.5_wp * ( c_interp_1 + c_interp_2 )
                ENDDO
 !
 !--            Then set the values along the matching grid-lines  
-               IF  ( MOD( ifl(l), igsr ) == 0 )  THEN
-!
-!--            First interpolate to the flux point using the 3rd-order WS scheme
-                  f_interp_1 = c31 * workarrc_t(nw-1,m,l) + c32 * workarrc_t(nw,m,l) + c33 * workarrc_t(nw+1,m,l)                  
-!                  f(k,jfl(m):jfu(m),ifl(l)) = workarrc_t(nw,m,l)
-                  f(k,jfl(m):jfu(m),ifl(l)) = f_interp_1
+               IF ( MOD( ifl(ip), igsr ) == 0 )  THEN
+!
+!--               First interpolate to the flux point using the 3rd-order WS scheme
+                  c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)            &
+                             + c33 * workarr_t(kpw+1,jp,ip)                  
+                  child_array(kc,jfl(jp):jfu(jp),ifl(ip)) = c_interp_1
                ENDIF
 
@@ -4304,11 +4344,11 @@
 !
 !--      Finally, set the values along the last matching grid-line  
-         IF  ( MOD( ifl(icrw), igsr ) == 0 )  THEN
-            DO  m = jcsw, jcnw
+         IF  ( MOD( ifl(iprw), igsr ) == 0 )  THEN
+            DO  jp = jpsw, jpnw
 !
 !--            First interpolate to the flux point using the 3rd-order WS scheme
-               f_interp_1 = c31 * workarrc_t(nw-1,m,icrw) + c32 * workarrc_t(nw,m,icrw) + c33 * workarrc_t(nw+1,m,icrw)
-!               f(k,jfl(m):jfu(m),ifl(icrw)) = workarrc_t(nw,m,icrw)
-               f(k,jfl(m):jfu(m),ifl(icrw)) = f_interp_1
+               c_interp_1 = c31 * workarr_t(kpw-1,jp,iprw) + c32 * workarr_t(kpw,jp,iprw)           &
+                          + c33 * workarr_t(kpw+1,jp,iprw)
+               child_array(kc,jfl(jp):jfu(jp),ifl(iprw)) = c_interp_1
             ENDDO
          ENDIF
@@ -4318,8 +4358,8 @@
 !--      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))
+         IF ( ifl(iprw) < nxr )  THEN
+            DO  jp = jpsw, jpnw
+               DO  ic = ifl(iprw) + 1, nxr
+                  child_array(kc,jfl(jp):jfu(jp),ic) = child_array(kc,jfl(jp):jfu(jp),ifl(iprw))
                ENDDO
             ENDDO
@@ -4328,23 +4368,23 @@
       ELSE IF  ( var == 'v' )  THEN
 
-         DO  l = iclw, icrw
-            DO  m = jcsw, jcnw-1
+         DO  ip = iplw, iprw
+            DO  jp = jpsw, jpnw-1
 !
 !--            First interpolate to the flux point using the 3rd-order WS scheme
-               f_interp_1 = c31 * workarrc_t(nw-1,m,l)   + c32 * workarrc_t(nw,m,l)   + c33 * workarrc_t(nw+1,m,l)
-               f_interp_2 = c31 * workarrc_t(nw-1,m+1,l) + c32 * workarrc_t(nw,m+1,l) + c33 * workarrc_t(nw+1,m+1,l)
+               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip)   + c32 * workarr_t(kpw,jp,ip)             &
+                          + c33 * workarr_t(kpw+1,jp,ip)
+               c_interp_2 = c31 * workarr_t(kpw-1,jp+1,ip) + c32 * workarr_t(kpw,jp+1,ip)           &
+                          + c33 * workarr_t(kpw+1,jp+1,ip)
 !
 !--            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) )
-                  f(k,j,ifl(l):ifu(l)) = 0.5_wp * ( f_interp_1 + f_interp_2 )
+               DO  jc = jfl(jp), jfl(jp+1)          
+                  child_array(kc,jc,ifl(ip):ifu(ip)) = 0.5_wp * ( c_interp_1 + c_interp_2 )
                ENDDO
 !
 !--            Then set the values along the matching grid-lines  
-               IF  ( MOD( jfl(m), jgsr ) == 0 )  THEN
-                  f_interp_1 = c31 * workarrc_t(nw-1,m,l) + c32 * workarrc_t(nw,m,l) + c33 * workarrc_t(nw+1,m,l)
-!                  f(k,jfl(m),ifl(l):ifu(l)) = workarrc_t(nw,m,l)
-                  f(k,jfl(m),ifl(l):ifu(l)) = f_interp_1
+               IF ( MOD( jfl(jp), jgsr ) == 0 )  THEN
+                  c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)            &
+                             + c33 * workarr_t(kpw+1,jp,ip)
+                  child_array(kc,jfl(jp),ifl(ip):ifu(ip)) = c_interp_1
                ENDIF
                
@@ -4354,11 +4394,11 @@
 !
 !--      Finally, set the values along the last matching grid-line
-         IF  ( MOD( jfl(jcnw), jgsr ) == 0 )  THEN
-            DO  l = iclw, icrw
+         IF ( MOD( jfl(jpnw), jgsr ) == 0 )  THEN
+            DO  ip = iplw, iprw
 !
 !--            First interpolate to the flux point using the 3rd-order WS scheme
-               f_interp_1 = c31 * workarrc_t(nw-1,jcnw,l) + c32 * workarrc_t(nw,jcnw,l) + c33 * workarrc_t(nw+1,jcnw,l)
-!               f(k,jfl(jcnw),ifl(l):ifu(l)) = workarrc_t(nw,jcnw,l)
-               f(k,jfl(jcnw),ifl(l):ifu(l)) = f_interp_1
+               c_interp_1 = c31 * workarr_t(kpw-1,jpnw,ip) + c32 * workarr_t(kpw,jpnw,ip)           &
+                          + c33 * workarr_t(kpw+1,jpnw,ip)
+               child_array(kc,jfl(jpnw),ifl(ip):ifu(ip)) = c_interp_1
             ENDDO
          ENDIF
@@ -4368,8 +4408,8 @@
 !--      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))
+         IF  ( jfl(jpnw) < nyn )  THEN
+            DO  ip = iplw, iprw
+               DO  jc = jfl(jpnw)+1, nyn
+                  child_array(kc,jc,ifl(ip):ifu(ip)) = child_array(kc,jfl(jpnw),ifl(ip):ifu(ip))
                ENDDO
             ENDDO
@@ -4378,13 +4418,13 @@
       ELSE  ! any scalar variable
 
-         DO  l = iclw, icrw
-            DO  m = jcsw, jcnw
+         DO  ip = iplw, iprw
+            DO  jp = jpsw, jpnw
 !
 !--            First interpolate to the flux point using the 3rd-order WS scheme
-               f_interp_1 = c31 * workarrc_t(nw-1,m,l) + c32 * workarrc_t(nw,m,l) + c33 * workarrc_t(nw+1,m,l)
-               DO  i = ifl(l), ifu(l)
-                  DO  j = jfl(m), jfu(m)
-!                     f(k,j,i) = workarrc_t(nw,m,l)
-                     f(k,j,i) = f_interp_1
+               c_interp_1 = c31 * workarr_t(kpw-1,jp,ip) + c32 * workarr_t(kpw,jp,ip)               &
+                          + c33 * workarr_t(kpw+1,jp,ip)
+               DO  ic = ifl(ip), ifu(ip)
+                  DO  jc = jfl(jp), jfu(jp)
+                     child_array(kc,jc,ic) = c_interp_1
                   ENDDO
                ENDDO
@@ -4397,5 +4437,5 @@
 !--   Just fill up the redundant second ghost-node layer in case of var == w.
       IF ( var == 'w' )  THEN
-         f(nzt+1,:,:) = f(nzt,:,:)
+         child_array(nzt+1,:,:) = child_array(nzt,:,:)
       ENDIF
 
@@ -4404,5 +4444,6 @@
 
 
-   SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu, ijkfc, var )
+   SUBROUTINE pmci_anterp_tophat( child_array, parent_array, kct, ifl, ifu, jfl, jfu, kfl, kfu,     &
+                                  ijkfc, var )
 !
 !--   Anterpolation of internal-node values to be used as the parent-domain
@@ -4413,55 +4454,67 @@
       IMPLICIT NONE
 
-      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  f         !< Child-grid array
-      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)  ::  fc        !< Parent-grid array
-      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
-      INTEGER(iwp), INTENT(IN) ::  kct                                             !< Top boundary index for anterpolation along z
-      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 end 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 end index of child cells belonging to certain parent cell - z direction
+      INTEGER(iwp), INTENT(IN) ::  kct  !< Top boundary index for anterpolation along z
+      
+      INTEGER(iwp), DIMENSION(0:pg%nz+1,jpsa:jpna,ipla:ipra), INTENT(IN) ::  ijkfc  !< number of child grid points contributing
+                                                                                    !< to a parent grid box
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(ipla:ipra), INTENT(IN) ::  ifu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - x direction
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - y direction
+      INTEGER(iwp), DIMENSION(jpsa:jpna), INTENT(IN) ::  jfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - y direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfl  !< Indicates start index of child cells belonging to certain
+                                                              !< parent cell - z direction
+      INTEGER(iwp), DIMENSION(0:pg%nz+1), INTENT(IN) ::  kfu  !< Indicates end index of child cells belonging to certain
+                                                              !< parent cell - z direction
+
+      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  child_array   !< Child-grid array
+      REAL(wp), DIMENSION(0:pg%nz+1,jps:jpn,ipl:ipr), INTENT(INOUT)  ::  parent_array  !< Parent-grid array
+
       CHARACTER(LEN=*), INTENT(IN) ::  var                   !< Variable symbol: 'u', 'v', 'w' or 's'
 !
 !--   Local variables:  
+      INTEGER(iwp) ::  ic              !< Running index x-direction - child grid
+      INTEGER(iwp) ::  ipl_anterp      !< Left boundary index for anterpolation along x
+      INTEGER(iwp) ::  ipr_anterp      !< Right boundary index for anterpolation along x
+      INTEGER(iwp) ::  jc              !< Running index y-direction - child grid
+      INTEGER(iwp) ::  jpn_anterp      !< North boundary index for anterpolation along y
+      INTEGER(iwp) ::  jps_anterp      !< South boundary index for anterpolation along y
+      INTEGER(iwp) ::  kc              !< Running index z-direction - child grid     
+      INTEGER(iwp) ::  kpb_anterp = 0  !< Bottom boundary index for anterpolation along z
+      INTEGER(iwp) ::  kpt_anterp      !< Top boundary index for anterpolation along z
+      INTEGER(iwp) ::  ip              !< Running index x-direction - parent grid
+      INTEGER(iwp) ::  jp              !< Running index y-direction - parent grid
+      INTEGER(iwp) ::  kp              !< Running index z-direction - parent grid
+      INTEGER(iwp) ::  var_flag        !< bit number used to flag topography on respective grid
+
       REAL(wp) ::  cellsum       !< sum of respective child cells belonging to parent cell 
-      INTEGER(iwp) ::  i         !< Running index x-direction - child grid
-      INTEGER(iwp) ::  iclant    !< Left boundary index for anterpolation along x
-      INTEGER(iwp) ::  icrant    !< Right boundary index for anterpolation along x
-      INTEGER(iwp) ::  j         !< Running index y-direction - child grid
-      INTEGER(iwp) ::  jcnant    !< North boundary index for anterpolation along y
-      INTEGER(iwp) ::  jcsant    !< South boundary index for anterpolation along y
-      INTEGER(iwp) ::  k         !< Running index z-direction - child grid     
-      INTEGER(iwp) ::  kcb = 0   !< Bottom boundary index for anterpolation along z
-      INTEGER(iwp) ::  kctant    !< Top boundary index for anterpolation along z
-      INTEGER(iwp) ::  l         !< Running index x-direction - parent grid
-      INTEGER(iwp) ::  m         !< Running index y-direction - parent grid
-      INTEGER(iwp) ::  n         !< Running index z-direction - parent grid
-      INTEGER(iwp) ::  var_flag  !< bit number used to flag topography on respective grid
-
-!
-!--   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_define_index_mapping.
+
+!
+!--   Define the index bounds ipl_anterp, ipr_anterp, jps_anterp and jpn_anterp.
+!--   Note that kcb_anterp is simply zero and kct_anterp enters here as a
+!--   parameter and it is determined in pmci_define_index_mapping.
 !--   Note that the 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. An additional buffer is
+!--   child) are always excluded from anterpolation, e.g. anterpolation is maximally
+!--   only from nzb:kct-1, as kct is used for interpolation. Similar restriction is
+!--   applied to the lateral boundaries as well. An additional buffer is
 !--   also applied (default value for anterpolation_buffer_width = 2) in order
 !--   to avoid unphysical accumulation of kinetic energy.
-      iclant = icl
-      icrant = icr
-      jcsant = jcs
-      jcnant = jcn
-!      kctant = kct - 1
-      kctant = kct - 1 - anterpolation_buffer_width
-      kcb  = 0
+      ipl_anterp = ipl
+      ipr_anterp = ipr
+      jps_anterp = jps
+      jpn_anterp = jpn
+      kpb_anterp = 0
+      kpt_anterp = kct - 1 - anterpolation_buffer_width
+
       IF ( nesting_mode /= 'vertical' )  THEN
 !
 !--      Set the anterpolation buffers on the lateral boundaries
-         iclant = MAX( icl, iplg + 3 + anterpolation_buffer_width )
-         icrant = MIN( icr, iprg - 3 - anterpolation_buffer_width )
-         jcsant = MAX( jcs, jpsg + 3 + anterpolation_buffer_width )
-         jcnant = MIN( jcn, jpng - 3 - anterpolation_buffer_width )
+         ipl_anterp = MAX( ipl, iplg + 3 + anterpolation_buffer_width )
+         ipr_anterp = MIN( ipr, iprg - 3 - anterpolation_buffer_width )
+         jps_anterp = MAX( jps, jpsg + 3 + anterpolation_buffer_width )
+         jpn_anterp = MIN( jpn, jpng - 3 - anterpolation_buffer_width )
          
       ENDIF
@@ -4478,18 +4531,18 @@
       ENDIF
 !
-!--   Note that l, m, and n are parent-grid indices and i,j, and k 
+!--   Note that ip, jp, and kp are parent-grid indices and ic,jc, and kc 
 !--   are child-grid indices.
-      DO  l = iclant, icrant
-         DO  m = jcsant, jcnant
+      DO  ip = ipl_anterp, ipr_anterp
+         DO  jp = jps_anterp, jpn_anterp
 !
 !--         For simplicity anterpolate within buildings and under elevated
 !--         terrain too
-            DO  n = kcb, kctant
+            DO  kp = kpb_anterp, kpt_anterp
                cellsum = 0.0_wp
-               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 ) )
+               DO  ic = ifl(ip), ifu(ip)
+                  DO  jc = jfl(jp), jfu(jp)
+                     DO  kc = kfl(kp), kfu(kp)
+                        cellsum = cellsum + MERGE( child_array(kc,jc,ic), 0.0_wp,                   &
+                             BTEST( wall_flags_0(kc,jc,ic), var_flag ) )
                      ENDDO
                   ENDDO
@@ -4501,7 +4554,6 @@
 !--            particular for the temperature. Therefore, in case cellsum is 
 !--            zero, keep the parent solution at this point. 
-
-               IF ( ijkfc(n,m,l) /= 0 )  THEN
-                  fc(n,m,l) = cellsum / REAL( ijkfc(n,m,l), KIND=wp )
+               IF ( ijkfc(kp,jp,ip) /= 0 )  THEN
+                  parent_array(kp,jp,ip) = cellsum / REAL( ijkfc(kp,jp,ip), KIND=wp )
                ENDIF
 
@@ -4539,13 +4591,13 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i   !< Index along x-direction
-    INTEGER(iwp) ::  ib  !< running index for aerosol size bins
-    INTEGER(iwp) ::  ic  !< running index for aerosol mass bins
-    INTEGER(iwp) ::  ig  !< running index for salsa gases
-    INTEGER(iwp) ::  j   !< Index along y-direction
-    INTEGER(iwp) ::  k   !< Index along z-direction
+    INTEGER(iwp) ::  ic  !< Index along x-direction
+    INTEGER(iwp) ::  jc  !< Index along y-direction
+    INTEGER(iwp) ::  kc  !< Index along z-direction
+    INTEGER(iwp) ::  lb  !< Running index for aerosol size bins
+    INTEGER(iwp) ::  lc  !< Running index for aerosol mass bins
+    INTEGER(iwp) ::  lg  !< Running index for salsa gases
     INTEGER(iwp) ::  m   !< Running index for surface type
-    INTEGER(iwp) ::  n   !< running index for number of chemical species
-
+    INTEGER(iwp) ::  n   !< Running index for number of chemical species
+    
 !
 !-- Set Dirichlet boundary conditions for horizontal velocity components
@@ -4554,18 +4606,18 @@
 !--    Upward-facing surfaces
        DO  m = 1, bc_h(0)%ns
-          i = bc_h(0)%i(m)            
-          j = bc_h(0)%j(m)
-          k = bc_h(0)%k(m)
-          u(k-1,j,i) = 0.0_wp
-          v(k-1,j,i) = 0.0_wp
+          ic = bc_h(0)%i(m)            
+          jc = bc_h(0)%j(m)
+          kc = bc_h(0)%k(m)
+          u(kc-1,jc,ic) = 0.0_wp
+          v(kc-1,jc,ic) = 0.0_wp
        ENDDO
 !
 !--    Downward-facing surfaces
        DO  m = 1, bc_h(1)%ns
-          i = bc_h(1)%i(m)            
-          j = bc_h(1)%j(m)
-          k = bc_h(1)%k(m)
-          u(k+1,j,i) = 0.0_wp
-          v(k+1,j,i) = 0.0_wp
+          ic = bc_h(1)%i(m)            
+          jc = bc_h(1)%j(m)
+          kc = bc_h(1)%k(m)
+          u(kc+1,jc,ic) = 0.0_wp
+          v(kc+1,jc,ic) = 0.0_wp
        ENDDO
     ENDIF
@@ -4574,16 +4626,16 @@
 !-- Upward-facing surfaces
     DO  m = 1, bc_h(0)%ns
-       i = bc_h(0)%i(m)            
-       j = bc_h(0)%j(m)
-       k = bc_h(0)%k(m)
-       w(k-1,j,i) = 0.0_wp
+       ic = bc_h(0)%i(m)            
+       jc = bc_h(0)%j(m)
+       kc = bc_h(0)%k(m)
+       w(kc-1,jc,ic) = 0.0_wp
     ENDDO
 !
 !-- Downward-facing surfaces
     DO  m = 1, bc_h(1)%ns
-       i = bc_h(1)%i(m)            
-       j = bc_h(1)%j(m)
-       k = bc_h(1)%k(m)
-       w(k+1,j,i) = 0.0_wp
+       ic = bc_h(1)%i(m)            
+       jc = bc_h(1)%j(m)
+       kc = bc_h(1)%k(m)
+       w(kc+1,jc,ic) = 0.0_wp
     ENDDO
 !
@@ -4592,14 +4644,14 @@
        IF ( ibc_pt_b == 1 )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             pt(k-1,j,i) = pt(k,j,i)
+             ic = bc_h(0)%i(m)            
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             pt(kc-1,jc,ic) = pt(kc,jc,ic)
           ENDDO
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)            
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-             pt(k+1,j,i) = pt(k,j,i)
+             ic = bc_h(1)%i(m)            
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+             pt(kc+1,jc,ic) = pt(kc,jc,ic)
           ENDDO   
        ENDIF
@@ -4610,31 +4662,31 @@
        IF ( ibc_q_b == 1 )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             q(k-1,j,i) = q(k,j,i)
+             ic = bc_h(0)%i(m)            
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             q(kc-1,jc,ic) = q(kc,jc,ic)
           ENDDO  
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)            
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-             q(k+1,j,i) = q(k,j,i)
+             ic = bc_h(1)%i(m)            
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+             q(kc+1,jc,ic) = q(kc,jc,ic)
           ENDDO  
        ENDIF
        IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             nc(k-1,j,i) = 0.0_wp
-             qc(k-1,j,i) = 0.0_wp
+             ic = bc_h(0)%i(m)            
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             nc(kc-1,jc,ic) = 0.0_wp
+             qc(kc-1,jc,ic) = 0.0_wp
           ENDDO  
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)            
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-
-             nc(k+1,j,i) = 0.0_wp
-             qc(k+1,j,i) = 0.0_wp
+             ic = bc_h(1)%i(m)            
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+
+             nc(kc+1,jc,ic) = 0.0_wp
+             qc(kc+1,jc,ic) = 0.0_wp
           ENDDO  
        ENDIF
@@ -4642,16 +4694,16 @@
        IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             nr(k-1,j,i) = 0.0_wp
-             qr(k-1,j,i) = 0.0_wp
+             ic = bc_h(0)%i(m)            
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             nr(kc-1,jc,ic) = 0.0_wp
+             qr(kc-1,jc,ic) = 0.0_wp
           ENDDO  
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)            
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-             nr(k+1,j,i) = 0.0_wp
-             qr(k+1,j,i) = 0.0_wp
+             ic = bc_h(1)%i(m)            
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+             nr(kc+1,jc,ic) = 0.0_wp
+             qr(kc+1,jc,ic) = 0.0_wp
           ENDDO  
        ENDIF
@@ -4663,14 +4715,14 @@
        IF ( ibc_s_b == 1 )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)            
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             s(k-1,j,i) = s(k,j,i)
+             ic = bc_h(0)%i(m)            
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             s(kc-1,jc,ic) = s(kc,jc,ic)
           ENDDO 
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)            
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-             s(k+1,j,i) = s(k,j,i)
+             ic = bc_h(1)%i(m)            
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+             s(kc+1,jc,ic) = s(kc,jc,ic)
           ENDDO  
        ENDIF
@@ -4682,14 +4734,14 @@
           DO  n = 1, nspec
              DO  m = 1, bc_h(0)%ns
-                i = bc_h(0)%i(m)            
-                j = bc_h(0)%j(m)
-                k = bc_h(0)%k(m)
-                chem_species(n)%conc(k-1,j,i) = chem_species(n)%conc(k,j,i)
+                ic = bc_h(0)%i(m)            
+                jc = bc_h(0)%j(m)
+                kc = bc_h(0)%k(m)
+                chem_species(n)%conc(kc-1,jc,ic) = chem_species(n)%conc(kc,jc,ic)
              ENDDO 
              DO  m = 1, bc_h(1)%ns
-                i = bc_h(1)%i(m)            
-                j = bc_h(1)%j(m)
-                k = bc_h(1)%k(m)
-                chem_species(n)%conc(k+1,j,i) = chem_species(n)%conc(k,j,i)
+                ic = bc_h(1)%i(m)            
+                jc = bc_h(1)%j(m)
+                kc = bc_h(1)%k(m)
+                chem_species(n)%conc(kc+1,jc,ic) = chem_species(n)%conc(kc,jc,ic)
              ENDDO
           ENDDO
@@ -4701,32 +4753,32 @@
        IF ( ibc_salsa_b == 1 )  THEN
           DO  m = 1, bc_h(0)%ns
-             i = bc_h(0)%i(m)
-             j = bc_h(0)%j(m)
-             k = bc_h(0)%k(m)
-             DO  ib = 1, nbins_aerosol
-                aerosol_number(ib)%conc(k-1,j,i) = aerosol_number(ib)%conc(k,j,i)
+             ic = bc_h(0)%i(m)
+             jc = bc_h(0)%j(m)
+             kc = bc_h(0)%k(m)
+             DO  lb = 1, nbins_aerosol
+                aerosol_number(lb)%conc(kc-1,jc,ic) = aerosol_number(lb)%conc(kc,jc,ic)
              ENDDO
-             DO  ic = 1, nbins_aerosol * ncomponents_mass
-                aerosol_mass(ic)%conc(k-1,j,i) = aerosol_mass(ic)%conc(k,j,i)
+             DO  lc = 1, nbins_aerosol * ncomponents_mass
+                aerosol_mass(lc)%conc(kc-1,jc,ic) = aerosol_mass(lc)%conc(kc,jc,ic)
              ENDDO
              IF ( .NOT. salsa_gases_from_chem )  THEN
-                DO  ig = 1, ngases_salsa
-                   salsa_gas(ig)%conc(k-1,j,i) = salsa_gas(ig)%conc(k,j,i)
+                DO  lg = 1, ngases_salsa
+                   salsa_gas(lg)%conc(kc-1,jc,ic) = salsa_gas(lg)%conc(kc,jc,ic)
                 ENDDO
              ENDIF
           ENDDO
           DO  m = 1, bc_h(1)%ns
-             i = bc_h(1)%i(m)
-             j = bc_h(1)%j(m)
-             k = bc_h(1)%k(m)
-             DO  ib = 1, nbins_aerosol
-                aerosol_number(ib)%conc(k+1,j,i) = aerosol_number(ib)%conc(k,j,i)
+             ic = bc_h(1)%i(m)
+             jc = bc_h(1)%j(m)
+             kc = bc_h(1)%k(m)
+             DO  lb = 1, nbins_aerosol
+                aerosol_number(lb)%conc(kc+1,jc,ic) = aerosol_number(lb)%conc(kc,jc,ic)
              ENDDO
-             DO  ic = 1, nbins_aerosol * ncomponents_mass
-                aerosol_mass(ic)%conc(k+1,j,i) = aerosol_mass(ic)%conc(k,j,i)
+             DO  lc = 1, nbins_aerosol * ncomponents_mass
+                aerosol_mass(lc)%conc(kc+1,jc,ic) = aerosol_mass(lc)%conc(kc,jc,ic)
              ENDDO
              IF ( .NOT. salsa_gases_from_chem )  THEN
-                DO  ig = 1, ngases_salsa
-                   salsa_gas(ig)%conc(k+1,j,i) = salsa_gas(ig)%conc(k,j,i)
+                DO  lg = 1, ngases_salsa
+                   salsa_gas(lg)%conc(kc+1,jc,ic) = salsa_gas(lg)%conc(kc,jc,ic)
                 ENDDO
              ENDIF
@@ -4736,5 +4788,4 @@
 
  END SUBROUTINE pmci_boundary_conds
-
-
+   
 END MODULE pmc_interface
