Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1789)
+++ palm/trunk/SOURCE/Makefile	(revision 1791)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# dependencies of header changed
 # 
 # Former revisions:
@@ -381,6 +381,6 @@
 global_min_max.o: modules.o mod_kinds.o
 header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model.o\
-   plant_canopy_model.o pmc_interface.o radiation_model.o spectrum.o \
-   subsidence.o 
+   plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \
+   radiation_model.o spectrum.o subsidence.o 
 impact_of_latent_heat.o: modules.o mod_kinds.o
 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1789)
+++ palm/trunk/SOURCE/header.f90	(revision 1791)
@@ -19,6 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
-! 
+! output of nesting informations of all domains
+!
 ! Former revisions:
 ! -----------------
@@ -284,7 +284,9 @@
                plant_canopy
 
+    USE pmc_handle_communicator,                                               &
+        ONLY:  pmc_get_model_info
+
     USE pmc_interface,                                                         &
-        ONLY:  cpl_id, cpl_parent_id, cpl_name, lower_left_coord_x,            &
-               lower_left_coord_y, nested_run, nesting_mode
+        ONLY:  nested_run, nesting_mode
 
     USE radiation_model_mod,                                                   &
@@ -311,4 +313,6 @@
     
     CHARACTER (LEN=26) ::  ver_rev             !<
+
+    CHARACTER (LEN=32) ::  cpl_name            !<
     
     CHARACTER (LEN=40) ::  output_format       !<
@@ -338,31 +342,38 @@
     CHARACTER (LEN=1), DIMENSION(1:3) ::  dir = (/ 'x', 'y', 'z' /)  !<
 
-    INTEGER(iwp) ::  av        !<
-    INTEGER(iwp) ::  bh        !<
-    INTEGER(iwp) ::  blx       !<
-    INTEGER(iwp) ::  bly       !<
-    INTEGER(iwp) ::  bxl       !<
-    INTEGER(iwp) ::  bxr       !<
-    INTEGER(iwp) ::  byn       !<
-    INTEGER(iwp) ::  bys       !<
-    INTEGER(iwp) ::  ch        !<
-    INTEGER(iwp) ::  count     !<
-    INTEGER(iwp) ::  cwx       !<
-    INTEGER(iwp) ::  cwy       !<
-    INTEGER(iwp) ::  cxl       !<
-    INTEGER(iwp) ::  cxr       !<
-    INTEGER(iwp) ::  cyn       !<
-    INTEGER(iwp) ::  cys       !<
-    INTEGER(iwp) ::  dim       !<
-    INTEGER(iwp) ::  i         !<
-    INTEGER(iwp) ::  io        !<
-    INTEGER(iwp) ::  j         !<
-    INTEGER(iwp) ::  k         !<
-    INTEGER(iwp) ::  l         !<
-    INTEGER(iwp) ::  ll        !<
-    INTEGER(iwp) ::  mpi_type  !<
+    INTEGER(iwp) ::  av             !<
+    INTEGER(iwp) ::  bh             !<
+    INTEGER(iwp) ::  blx            !<
+    INTEGER(iwp) ::  bly            !<
+    INTEGER(iwp) ::  bxl            !<
+    INTEGER(iwp) ::  bxr            !<
+    INTEGER(iwp) ::  byn            !<
+    INTEGER(iwp) ::  bys            !<
+    INTEGER(iwp) ::  ch             !<
+    INTEGER(iwp) ::  count          !<
+    INTEGER(iwp) ::  cpl_parent_id  !<
+    INTEGER(iwp) ::  cwx            !<
+    INTEGER(iwp) ::  cwy            !<
+    INTEGER(iwp) ::  cxl            !<
+    INTEGER(iwp) ::  cxr            !<
+    INTEGER(iwp) ::  cyn            !<
+    INTEGER(iwp) ::  cys            !<
+    INTEGER(iwp) ::  dim            !<
+    INTEGER(iwp) ::  i              !<
+    INTEGER(iwp) ::  io             !<
+    INTEGER(iwp) ::  j              !<
+    INTEGER(iwp) ::  k              !<
+    INTEGER(iwp) ::  l              !<
+    INTEGER(iwp) ::  ll             !<
+    INTEGER(iwp) ::  mpi_type       !<
+    INTEGER(iwp) ::  my_cpl_id      !<
+    INTEGER(iwp) ::  n              !<
+    INTEGER(iwp) ::  ncpl           !<
+    INTEGER(iwp) ::  npe_total      !<
     
     REAL(wp) ::  canopy_height                    !< canopy height (in m)
     REAL(wp) ::  cpuseconds_per_simulated_second  !<
+    REAL(wp) ::  lower_left_coord_x               !< x-coordinate of nest domain
+    REAL(wp) ::  lower_left_coord_y               !< y-coordinate of nest domain
 
 !
@@ -468,6 +479,23 @@
 !-- Nesting informations
     IF ( nested_run )  THEN
-       WRITE ( io, 600 )  cpl_id, TRIM( cpl_name ), cpl_parent_id,             &
-                          nesting_mode, lower_left_coord_x, lower_left_coord_y
+
+       WRITE ( io, 600 )  TRIM( nesting_mode )
+       CALL pmc_get_model_info( ncpl = ncpl, cpl_id = my_cpl_id )
+
+       DO  n = 1, ncpl
+          CALL pmc_get_model_info( request_for_cpl_id = n, cpl_name = cpl_name,&
+                                   cpl_parent_id = cpl_parent_id,              &
+                                   lower_left_x = lower_left_coord_x,          &
+                                   lower_left_y = lower_left_coord_y,          &
+                                   npe_total = npe_total )
+          IF ( n == my_cpl_id )  THEN
+             char1 = '*'
+          ELSE
+             char1 = ' '
+          ENDIF
+          WRITE ( io, 601 )  TRIM( char1 ), n, cpl_parent_id, npe_total,       &
+                             lower_left_coord_x, lower_left_coord_y,           &
+                             TRIM( cpl_name )
+       ENDDO
     ENDIF
     WRITE ( io, 99 )
@@ -2443,8 +2471,10 @@
 513 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order ' // & 
             '+ monotonic adjustment')
-600 FORMAT (/' Nesting informations:'/                                        &
-            ' Nest id / name:                   ',I2.2,' / ',A,' (parent id: ',I2.2,')'/ &
-            ' Nesting mode:                     ',A/ &
-            ' Lower left corner coordinates:    ','x = ',F8.2,' m, y = ',F8.2,' m'/)
+600 FORMAT (/' Nesting informations:'/ &
+            ' --------------------'/ &
+            ' Nesting mode:                     ',A// &
+            ' Nest id  parent  number   lower left coordinates   name'/ &
+            ' (*=me)     id    of PEs      x (m)     y (m)' )
+601 FORMAT (2X,A1,1X,I2.2,6X,I2.2,5X,I5,5X,F8.2,2X,F8.2,5X,A)
 
  END SUBROUTINE header
Index: palm/trunk/SOURCE/pmc_client.f90
===================================================================
--- palm/trunk/SOURCE/pmc_client.f90	(revision 1789)
+++ palm/trunk/SOURCE/pmc_client.f90	(revision 1791)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Debug write-statement commented out
 !
 ! Former revisions:
@@ -150,5 +150,5 @@
         end do
 
-        if(me%model_rank == 0) write(0,'(a,5i6)') 'PMC_ClientInit ',me%model_rank,me%model_npes,me%inter_npes,me%intra_rank
+!        if(me%model_rank == 0) write(0,'(a,5i6)') 'PMC_ClientInit ',me%model_rank,me%model_npes,me%inter_npes,me%intra_rank
 
         return
Index: palm/trunk/SOURCE/pmc_handle_communicator.f90
===================================================================
--- palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1789)
+++ palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1791)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-!
+! m_nrofcpl renamed m_ncpl,
+! pmc_get_local_model_info renamed pmc_get_model_info, some keywords also
+! renamed and some added,
+! debug write-statements commented out
 !
 ! Former revisions:
@@ -86,6 +89,7 @@
    INTEGER                                    :: m_my_CPL_id  !Coupler id of this model
    INTEGER                                    :: m_Parent_id  !Coupler id of parent of this model
-   INTEGER                                    :: m_NrOfCpl    !Number of Coupler in layout file
-   TYPE(PMC_layout),DIMENSION(PMC_MAX_MODELL) :: m_couplers   !Information of all coupler
+   INTEGER                                    :: m_ncpl       !Number of Couplers in layout file
+
+   TYPE(PMC_layout),DIMENSION(PMC_MAX_MODELL) :: m_couplers   !Information of all couplers
 
    ! MPI settings
@@ -110,9 +114,9 @@
    END INTERFACE pmc_is_rootmodel
 
-   INTERFACE PMC_get_local_model_info
-      MODULE PROCEDURE PMC_get_local_model_info
-   END INTERFACE PMC_get_local_model_info
-
-   PUBLIC pmc_get_local_model_info, pmc_init_model, pmc_is_rootmodel
+   INTERFACE pmc_get_model_info
+      MODULE PROCEDURE pmc_get_model_info
+   END INTERFACE pmc_get_model_info
+
+   PUBLIC pmc_get_model_info, pmc_init_model, pmc_is_rootmodel
 
  CONTAINS
@@ -160,5 +164,5 @@
 !--         Calculate start PE of every model
             start_pe(1) = 0
-            DO  i = 2, m_nrofcpl+1
+            DO  i = 2, m_ncpl+1
                start_pe(i) = start_pe(i-1) + m_couplers(i-1)%npe_total
             ENDDO
@@ -167,8 +171,8 @@
 !--         The number of cores provided with the run must be the same as the
 !--         total sum of cores required by all nest domains
-            IF ( start_pe(m_nrofcpl+1) /= m_world_npes )  THEN
+            IF ( start_pe(m_ncpl+1) /= m_world_npes )  THEN
                WRITE ( message_string, '(A,I6,A,I6,A)' )                       &
                                'nesting-setup requires more MPI procs (',      &
-                               start_pe(m_nrofcpl+1), ') than provided (',     &
+                               start_pe(m_ncpl+1), ') than provided (',        &
                                m_world_npes,')'
                CALL message( 'pmc_init_model', 'PA0229', 3, 2, 0, 6, 0 )
@@ -202,10 +206,10 @@
       ENDIF
 
-      CALL MPI_BCAST( m_nrofcpl, 1,          MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
-      CALL MPI_BCAST( start_pe, m_nrofcpl+1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
+      CALL MPI_BCAST( m_ncpl,          1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
+      CALL MPI_BCAST( start_pe, m_ncpl+1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
 
 !
 !--   Broadcast coupling layout
-      DO  i = 1, m_nrofcpl
+      DO  i = 1, m_ncpl
          CALL MPI_BCAST( m_couplers(i)%name, LEN( m_couplers(i)%name ), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
          CALL MPI_BCAST( m_couplers(i)%id,           1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat )
@@ -219,5 +223,5 @@
 !
 !--   Assign global MPI processes to individual models by setting the couple id
-      DO  i = 1, m_nrofcpl
+      DO  i = 1, m_ncpl
          IF ( m_world_rank >= start_pe(i)  .AND.  m_world_rank < start_pe(i+1) ) &
          THEN
@@ -241,5 +245,5 @@
 !
 !--   Broadcast (from PE 0) the parent id and id of every model
-      DO  i = 1, m_nrofcpl
+      DO  i = 1, m_ncpl
          CALL MPI_BCAST( m_couplers(i)%parent_id, 1, MPI_INTEGER, 0,           &
                          MPI_COMM_WORLD, istat )
@@ -257,5 +261,5 @@
 !--   different colors.
 !--   The grouping was done above with MPI_COMM_SPLIT
-      DO  i = 2, m_nrofcpl
+      DO  i = 2, m_ncpl
 
          IF ( m_couplers(i)%parent_id == m_my_cpl_id )  THEN
@@ -292,5 +296,5 @@
 
       clientcount = 0
-      DO  i = 2, m_nrofcpl
+      DO  i = 2, m_ncpl
          IF ( activeserver(i) == 1 )  THEN
             clientcount = clientcount + 1
@@ -321,7 +325,8 @@
 
 !
-!-- Make module private variables available to palm
-   SUBROUTINE pmc_get_local_model_info( my_cpl_id, my_cpl_parent_id, cpl_name, &
-                                        npe_total, lower_left_x, lower_left_y )
+!-- Provide module private variables of the pmc for PALM
+    SUBROUTINE pmc_get_model_info( cpl_id, cpl_name, cpl_parent_id,            &
+                                   lower_left_x, lower_left_y, ncpl, npe_total,&
+                                   request_for_cpl_id )
 
       USE kinds
@@ -330,18 +335,53 @@
 
       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL ::  cpl_name
-      INTEGER, INTENT(OUT), OPTIONAL          ::  my_cpl_id
-      INTEGER, INTENT(OUT), OPTIONAL          ::  my_cpl_parent_id
-      INTEGER, INTENT(OUT), OPTIONAL          ::  npe_total
-      REAL(wp), INTENT(OUT), OPTIONAL         ::  lower_left_x
-      REAL(wp), INTENT(OUT), OPTIONAL         ::  lower_left_y
-
-      IF ( PRESENT( my_cpl_id )           )  my_cpl_id        = m_my_cpl_id
-      IF ( PRESENT( my_cpl_parent_id )    )  my_cpl_parent_id = m_couplers(my_cpl_id)%parent_id
-      IF ( PRESENT( cpl_name )            )  cpl_name         = m_couplers(my_cpl_id)%name
-      IF ( PRESENT( npe_total )           )  npe_total        = m_couplers(my_cpl_id)%npe_total
-      IF ( PRESENT( lower_left_x )        )  lower_left_x     = m_couplers(my_cpl_id)%lower_left_x
-      IF ( PRESENT( lower_left_y )        )  lower_left_y     = m_couplers(my_cpl_id)%lower_left_y
-
-   END SUBROUTINE pmc_get_local_model_info
+
+      INTEGER, INTENT(IN), OPTIONAL ::  request_for_cpl_id
+
+      INTEGER, INTENT(OUT), OPTIONAL ::  cpl_id
+      INTEGER, INTENT(OUT), OPTIONAL ::  cpl_parent_id
+      INTEGER, INTENT(OUT), OPTIONAL ::  ncpl
+      INTEGER, INTENT(OUT), OPTIONAL ::  npe_total
+
+      INTEGER ::  requested_cpl_id
+
+      REAL(wp), INTENT(OUT), OPTIONAL ::  lower_left_x
+      REAL(wp), INTENT(OUT), OPTIONAL ::  lower_left_y
+
+!
+!--   Set the requested coupler id
+      IF ( PRESENT( request_for_cpl_id ) )  THEN
+         requested_cpl_id = request_for_cpl_id
+!
+!--      Check for allowed range of values
+         IF ( requested_cpl_id < 1 .OR. requested_cpl_id > m_ncpl )  RETURN
+      ELSE
+         requested_cpl_id = m_my_cpl_id
+      ENDIF
+
+!
+!--   Return the requested information
+      IF ( PRESENT( cpl_id )        )  THEN
+         cpl_id        = requested_cpl_id
+      ENDIF
+      IF ( PRESENT( cpl_parent_id ) )  THEN
+         cpl_parent_id = m_couplers(requested_cpl_id)%parent_id
+      ENDIF
+      IF ( PRESENT( cpl_name )      )  THEN
+         cpl_name      = m_couplers(requested_cpl_id)%name
+      ENDIF
+      IF ( PRESENT( ncpl )          )  THEN
+         ncpl          = m_ncpl
+      ENDIF
+      IF ( PRESENT( npe_total )     )  THEN
+         npe_total     = m_couplers(requested_cpl_id)%npe_total
+      ENDIF
+      IF ( PRESENT( lower_left_x )  )  THEN
+         lower_left_x  = m_couplers(requested_cpl_id)%lower_left_x
+      ENDIF
+      IF ( PRESENT( lower_left_y )  )  THEN
+         lower_left_y  = m_couplers(requested_cpl_id)%lower_left_y
+      ENDIF
+
+   END SUBROUTINE pmc_get_model_info
 
 
@@ -364,5 +404,5 @@
 
     INTEGER, INTENT(INOUT) ::  pmc_status
-    INTEGER                ::  i, istat, iunit
+    INTEGER                ::  i, istat
 
     TYPE(pmc_layout), DIMENSION(pmc_max_modell) ::  domain_layouts
@@ -374,6 +414,5 @@
 !-- Initialize some coupling variables
     domain_layouts(1:pmc_max_modell)%id = -1
-    m_nrofcpl =   0
-    iunit     = 345
+    m_ncpl =   0
 
     pmc_status = pmc_status_ok
@@ -412,17 +451,13 @@
 !-- Get the number of nested models given in the nestpar-NAMELIST
     DO  i = 1, pmc_max_modell
-
-       IF ( m_couplers(i)%id /= -1  .AND.  i <= pmc_max_modell )  THEN
-          WRITE ( 0, '(A,A,1X,3I7,1X,2F10.2)' )  'Set up Model  ',             &
-                             TRIM( m_couplers(i)%name ), m_couplers(i)%id,     &
-                             m_couplers(i)%Parent_id, m_couplers(i)%npe_total, &
-                             m_couplers(i)%lower_left_x,                       &
-                             m_couplers(i)%lower_left_y
-       ELSE
-!
-!--       When id=-1 is found for the first time, the list of domains is
-!--       finished (or latest after pmc_max_modell entries
-          m_nrofcpl = i - 1
-          EXIT
+!
+!--    When id=-1 is found for the first time, the list of domains is finished
+       IF ( m_couplers(i)%id == -1  .OR.  i == pmc_max_modell )  THEN
+          IF ( m_couplers(i)%id == -1 )  THEN
+             m_ncpl = i - 1
+             EXIT
+          ELSE
+             m_ncpl = pmc_max_modell
+          ENDIF
        ENDIF
 
Index: palm/trunk/SOURCE/pmc_interface.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface.f90	(revision 1789)
+++ palm/trunk/SOURCE/pmc_interface.f90	(revision 1791)
@@ -20,5 +20,10 @@
 ! Current revisions:
 ! ------------------
-!
+! routine pmci_update_new removed,
+! pmc_get_local_model_info renamed pmc_get_model_info, some keywords also
+! renamed,
+! filling up redundant ghost points introduced,
+! some index bound variables renamed,
+! further formatting cleanup
 !
 ! Former revisions:
@@ -111,5 +116,5 @@
 
     USE pmc_handle_communicator,                                               &
-        ONLY:  pmc_get_local_model_info, pmc_init_model, pmc_is_rootmodel,     &
+        ONLY:  pmc_get_model_info, pmc_init_model, pmc_is_rootmodel,           &
                pmc_no_namelist_found, pmc_server_for_client
 
@@ -128,155 +133,151 @@
     IMPLICIT NONE
 
-!-- TO_DO: a lot of lines (including comments) in this file exceed the 80 char
-!--        limit. Try to reduce as much as possible
-
-!-- TO_DO: are all of these variables following now really PUBLIC?
-!--        Klaus and I guess they are not
-    PRIVATE    !:  Note that the default publicity is here set to private.
+    PRIVATE
 
 !
 !-- Constants
-    INTEGER(iwp), PARAMETER, PUBLIC ::  client_to_server = 2   !:
-    INTEGER(iwp), PARAMETER, PUBLIC ::  server_to_client = 1   !:
+    INTEGER(iwp), PARAMETER ::  client_to_server = 2   !:
+    INTEGER(iwp), PARAMETER ::  server_to_client = 1   !:
 
 !
 !-- Coupler setup
-    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_id  = 1            !:
-    CHARACTER(LEN=32), PUBLIC, SAVE ::  cpl_name               !:
-    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_npe_total          !:
-    INTEGER(iwp), PUBLIC, SAVE      ::  cpl_parent_id          !:
+    INTEGER(iwp), SAVE      ::  cpl_id  = 1            !:
+    CHARACTER(LEN=32), SAVE ::  cpl_name               !:
+    INTEGER(iwp), SAVE      ::  cpl_npe_total          !:
+    INTEGER(iwp), SAVE      ::  cpl_parent_id          !:
 
 !
 !-- Control parameters, will be made input parameters later
-    CHARACTER(LEN=7), PUBLIC, SAVE ::  nesting_mode = 'two-way'  !: steering parameter for one- or two-way nesting
-
-    LOGICAL, PUBLIC, SAVE ::  nested_run = .FALSE.  !: general switch if nested run or not
-
-    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_l = -1.0_wp   !:
-    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_r = -1.0_wp   !:
-    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_s = -1.0_wp   !:
-    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_n = -1.0_wp   !:
-    REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_t = -1.0_wp   !:
+    CHARACTER(LEN=7), SAVE ::  nesting_mode = 'two-way'  !: steering parameter
+                                                         !: for 1- or 2-way nesting
+
+    LOGICAL, SAVE ::  nested_run = .FALSE.  !: general switch
+
+    REAL(wp), SAVE ::  anterp_relax_length_l = -1.0_wp   !:
+    REAL(wp), SAVE ::  anterp_relax_length_r = -1.0_wp   !:
+    REAL(wp), SAVE ::  anterp_relax_length_s = -1.0_wp   !:
+    REAL(wp), SAVE ::  anterp_relax_length_n = -1.0_wp   !:
+    REAL(wp), SAVE ::  anterp_relax_length_t = -1.0_wp   !:
 
 !
 !-- Geometry
-    REAL(wp), PUBLIC, SAVE                    ::  area_t               !:
+    REAL(wp), SAVE                            ::  area_t               !:
     REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_x              !:
     REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_y              !:
-    REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_x   !:
-    REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_y   !:
+    REAL(wp), SAVE                            ::  lower_left_coord_x   !:
+    REAL(wp), SAVE                            ::  lower_left_coord_y   !:
 
 !
 !-- Client coarse data arrays
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ec   !:
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ptc  !:
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  uc   !:
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  vc   !:
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  wc   !:
-    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  qc   !:
-
-    INTEGER(iwp), DIMENSION(5)                          ::  coarse_bound   !:
-    REAL(wp), PUBLIC, SAVE                              ::  xexl           !:
-    REAL(wp), PUBLIC, SAVE                              ::  xexr           !:
-    REAL(wp), PUBLIC, SAVE                              ::  yexs           !:
-    REAL(wp), PUBLIC, SAVE                              ::  yexn           !:
-    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !:
-    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !:
-    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !:
-    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !:
-    REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !:
-
-!
-!-- Client interpolation coefficients and client-array indices to be precomputed and stored.
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
-
-!
-!-- Client index arrays and log-ratio arrays for the log-law near-wall corrections.
-!-- These are not truly 3-D arrays but multiply 2-D arrays.
-    INTEGER(iwp), PUBLIC, SAVE :: ncorr  !: ncorr is the 4th dimension of the log_ratio-arrays
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
+    INTEGER(iwp), DIMENSION(5)                  ::  coarse_bound   !:
+
+    REAL(wp), SAVE                              ::  xexl           !:
+    REAL(wp), SAVE                              ::  xexr           !:
+    REAL(wp), SAVE                              ::  yexs           !:
+    REAL(wp), SAVE                              ::  yexn           !:
+    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !:
+    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !:
+    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !:
+    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !:
+    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !:
+
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec   !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc  !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uc   !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vc   !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc   !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc   !:
+
+!
+!-- Client interpolation coefficients and client-array indices to be precomputed
+!-- and stored.
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
+
+!
+!-- Client index arrays and log-ratio arrays for the log-law near-wall
+!-- corrections. These are not truly 3-D arrays but multiple 2-D arrays.
+    INTEGER(iwp), SAVE :: ncorr  !: 4th dimension of the log_ratio-arrays
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
 
 !
 !-- Upper bounds for k in anterpolation.
-    INTEGER(iwp), PUBLIC, SAVE ::  kceu   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  kcew   !:
+    INTEGER(iwp), SAVE ::  kctu   !:
+    INTEGER(iwp), SAVE ::  kctw   !:
 
 !
 !-- Upper bound for k in log-law correction in interpolation.
-    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_l   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_n   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_r   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_s   !:
+    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_l   !:
+    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_n   !:
+    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_r   !:
+    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_s   !:
 
 !
 !-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
-    INTEGER(iwp), PUBLIC, SAVE ::  nhll   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nhlr   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nhls   !:
-    INTEGER(iwp), PUBLIC, SAVE ::  nhln   !:
+    INTEGER(iwp), SAVE ::  nhll   !:
+    INTEGER(iwp), SAVE ::  nhlr   !:
+    INTEGER(iwp), SAVE ::  nhls   !:
+    INTEGER(iwp), SAVE ::  nhln   !:
 
 !
 !-- Spatial under-relaxation coefficients for anterpolation.
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
-    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
+    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
 
 !
 !-- Client-array indices to be precomputed and stored for anterpolation.
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
-    INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
-
-!
-!-- Module private variables.
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
+    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
+
     INTEGER(iwp), DIMENSION(3)          ::  define_coarse_grid_int    !:
     REAL(wp), DIMENSION(7)              ::  define_coarse_grid_real   !:
@@ -301,5 +302,5 @@
     END TYPE coarsegrid_def
                                           
-    TYPE(coarsegrid_def), SAVE             ::  cg   !:
+    TYPE(coarsegrid_def), SAVE ::  cg   !:
 
 
@@ -340,8 +341,8 @@
     END INTERFACE pmci_set_swaplevel
 
-    INTERFACE pmci_update_new
-       MODULE PROCEDURE pmci_update_new
-    END INTERFACE
-
+    PUBLIC anterp_relax_length_l, anterp_relax_length_r,                       &
+           anterp_relax_length_s, anterp_relax_length_n,                       &
+           anterp_relax_length_t, client_to_server, cpl_id, nested_run,        &
+           nesting_mode, server_to_client
     PUBLIC pmci_client_datatrans
     PUBLIC pmci_client_initialize
@@ -354,5 +355,4 @@
     PUBLIC pmci_server_synchronize
     PUBLIC pmci_set_swaplevel
-    PUBLIC pmci_update_new
 
 
@@ -393,13 +393,11 @@
 !-- Get some variables required by the pmc-interface (and in some cases in the
 !-- PALM code out of the pmci) out of the pmc-core
-    CALL pmc_get_local_model_info( my_cpl_id = cpl_id,                         &
-                                   my_cpl_parent_id = cpl_parent_id,           &
-                                   cpl_name = cpl_name,                        &
-                                   npe_total = cpl_npe_total,                  &
-                                   lower_left_x = lower_left_coord_x,          &
-                                   lower_left_y = lower_left_coord_y )
+    CALL pmc_get_model_info( cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,   &
+                             cpl_name = cpl_name, npe_total = cpl_npe_total,   &
+                             lower_left_x = lower_left_coord_x,                &
+                             lower_left_y = lower_left_coord_y )
 !
 !-- Set the steering switch which tells the models that they are nested (of
-!-- course the root domain (cpl_id = 1 ) is not nested)
+!-- course the root domain (cpl_id = 1) is not nested)
     IF ( cpl_id >= 2 )  THEN
        nest_domain = .TRUE.
@@ -435,7 +433,14 @@
 
     CALL location_message( 'setup the nested model configuration', .FALSE. )
-    CALL pmci_setup_coordinates   !:  Compute absolute coordinates valid for all models
-    CALL pmci_setup_client        !:  Initialize PMC Client (Must be called before pmc_setup_server)
-    CALL pmci_setup_server        !:  Initialize PMC Server
+!
+!-- Compute absolute coordinates for all models
+    CALL pmci_setup_coordinates
+!
+!-- Initialize the client (must be called before pmc_setup_server)
+    CALL pmci_setup_client
+!
+!-- Initialize PMC Server
+    CALL pmci_setup_server
+
     CALL location_message( 'finished', .TRUE. )
 
@@ -449,22 +454,26 @@
     IMPLICIT NONE
 
-    INTEGER(iwp)               ::  client_id        !:
-    INTEGER(iwp)               ::  ierr             !:
-    INTEGER(iwp)               ::  i                !:
-    INTEGER(iwp)               ::  j                !:
-    INTEGER(iwp)               ::  k                !:
-    INTEGER(iwp)               ::  m                !:
-    INTEGER(iwp)               ::  nomatch          !:
-    INTEGER(iwp)               ::  nx_cl            !:
-    INTEGER(iwp)               ::  ny_cl            !:
-    INTEGER(iwp)               ::  nz_cl            !:
-    INTEGER(iwp), DIMENSION(5) ::  val              !:
-    REAL(wp), DIMENSION(1)     ::  fval             !:
-    REAL(wp)                   ::  dx_cl            !:
-    REAL(wp)                   ::  dy_cl            !:
-    REAL(wp)                   ::  xez              !:
-    REAL(wp)                   ::  yez              !:
-    CHARACTER(len=32)          ::  mychannel
-    CHARACTER(len=32)          ::  myname
+    CHARACTER(LEN=32) ::  myname
+
+    INTEGER(iwp) ::  client_id        !:
+    INTEGER(iwp) ::  ierr             !:
+    INTEGER(iwp) ::  i                !:
+    INTEGER(iwp) ::  j                !:
+    INTEGER(iwp) ::  k                !:
+    INTEGER(iwp) ::  m                !:
+    INTEGER(iwp) ::  nomatch          !:
+    INTEGER(iwp) ::  nx_cl            !:
+    INTEGER(iwp) ::  ny_cl            !:
+    INTEGER(iwp) ::  nz_cl            !:
+
+    INTEGER(iwp), DIMENSION(5) ::  val    !:
+
+    REAL(wp) ::  dx_cl            !:
+    REAL(wp) ::  dy_cl            !:
+    REAL(wp) ::  xez              !:
+    REAL(wp) ::  yez              !:
+
+    REAL(wp), DIMENSION(1) ::  fval             !:
+
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !:
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !:
@@ -472,5 +481,5 @@
 
 !
-!   Initialize the PMC server
+!   Initialize the pmc server
     CALL pmc_serverinit
 
@@ -478,4 +487,5 @@
 !-- Get coordinates from all clients
     DO  m = 1, SIZE( pmc_server_for_client ) - 1
+
        client_id = pmc_server_for_client(m)
        IF ( myid == 0 )  THEN        
@@ -510,5 +520,5 @@
           CALL pmc_recv_from_client( client_id, cl_coord_y, SIZE( cl_coord_y ),&
                                      0, 12, ierr )
-          WRITE ( 0, * )  'receive from pmc Client ', client_id, nx_cl, ny_cl
+!          WRITE ( 0, * )  'receive from pmc Client ', client_id, nx_cl, ny_cl
          
           define_coarse_grid_real(1) = lower_left_coord_x
@@ -529,8 +539,11 @@
           xez = ( nbgp + 1 ) * dx 
           yez = ( nbgp + 1 ) * dy 
-          IF ( cl_coord_x(0) < define_coarse_grid_real(1) + xez )          nomatch = 1 
-          IF ( cl_coord_x(nx_cl + 1) > define_coarse_grid_real(5) - xez )  nomatch = 1
-          IF ( cl_coord_y(0) < define_coarse_grid_real(2) + yez )          nomatch = 1 
-          IF ( cl_coord_y(ny_cl + 1) > define_coarse_grid_real(6) - yez )  nomatch = 1
+          IF ( ( cl_coord_x(0) < define_coarse_grid_real(1) + xez )       .OR. &
+               ( cl_coord_x(nx_cl+1) > define_coarse_grid_real(5) - xez ) .OR. &
+               ( cl_coord_y(0) < define_coarse_grid_real(2) + yez )       .OR. &
+               ( cl_coord_y(ny_cl+1) > define_coarse_grid_real(6) - yez ) )    &
+          THEN
+             nomatch = 1
+          ENDIF
 
           DEALLOCATE( cl_coord_x )
@@ -539,21 +552,23 @@
 !
 !--       Send coarse grid information to client
-          CALL pmc_send_to_client( client_id, Define_coarse_grid_real,         &
-                                   SIZE(define_coarse_grid_real), 0,           &
-                                   21, ierr )
-          CALL pmc_send_to_client( client_id, Define_coarse_grid_int,  3, 0,   &
+          CALL pmc_send_to_client( client_id, define_coarse_grid_real,         &
+                                   SIZE( define_coarse_grid_real ), 0, 21,     &
+                                   ierr )
+          CALL pmc_send_to_client( client_id, define_coarse_grid_int,  3, 0,   &
                                    22, ierr )
 
 !
 !--       Send local grid to client
-          CALL pmc_send_to_client( client_id, coord_x, nx+1+2*nbgp, 0, 24, ierr )
-          CALL pmc_send_to_client( client_id, coord_y, ny+1+2*nbgp, 0, 25, ierr )
+          CALL pmc_send_to_client( client_id, coord_x, nx+1+2*nbgp, 0, 24,     &
+                                   ierr )
+          CALL pmc_send_to_client( client_id, coord_y, ny+1+2*nbgp, 0, 25,     &
+                                   ierr )
 
 !
 !--       Also send the dzu-, dzw-, zu- and zw-arrays here
-          CALL pmc_send_to_client( client_id, dzu, nz_cl + 1, 0, 26, ierr )
-          CALL pmc_send_to_client( client_id, dzw, nz_cl + 1, 0, 27, ierr )
-          CALL pmc_send_to_client( client_id, zu,  nz_cl + 2, 0, 28, ierr )
-          CALL pmc_send_to_client( client_id, zw,  nz_cl + 2, 0, 29, ierr )
+          CALL pmc_send_to_client( client_id, dzu, nz_cl+1, 0, 26, ierr )
+          CALL pmc_send_to_client( client_id, dzw, nz_cl+1, 0, 27, ierr )
+          CALL pmc_send_to_client( client_id, zu,  nz_cl+2, 0, 28, ierr )
+          CALL pmc_send_to_client( client_id, zw,  nz_cl+2, 0, 29, ierr )
 
        ENDIF
@@ -567,9 +582,12 @@
       
        CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
-       
+
+!
+!--    TO_DO: Klaus: please give a comment what is done here
        CALL pmci_create_index_list
 
 !
 !--    Include couple arrays into server content
+!--    TO_DO: Klaus: please give a more meaningful comment
        CALL pmc_s_clear_next_array_list
        DO  WHILE ( pmc_s_getnextarray( client_id, myname ) )
@@ -587,25 +605,27 @@
        IMPLICIT NONE
 
+       INTEGER(iwp) ::  i                  !:
+       INTEGER(iwp) ::  ic                 !:
+       INTEGER(iwp) ::  ierr               !:
+       INTEGER(iwp) ::  j                  !:
+       INTEGER(iwp) ::  k                  !:
+       INTEGER(iwp) ::  m                  !:
+       INTEGER(iwp) ::  n                  !:
+       INTEGER(iwp) ::  npx                !:
+       INTEGER(iwp) ::  npy                !:
+       INTEGER(iwp) ::  nrx                !:
+       INTEGER(iwp) ::  nry                !:
+       INTEGER(iwp) ::  px                 !:
+       INTEGER(iwp) ::  py                 !:
+       INTEGER(iwp) ::  server_pe          !:
+
+       INTEGER(iwp), DIMENSION(2) ::  scoord             !:
+       INTEGER(iwp), DIMENSION(2) ::  size_of_array      !:
+
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  coarse_bound_all   !:
-       INTEGER(iwp)                               ::  i                  !:
-       INTEGER(iwp)                               ::  ic                 !:
-       INTEGER(iwp)                               ::  ierr               !:
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list         !:
-       INTEGER(iwp)                               ::  j                  !:
-       INTEGER(iwp)                               ::  k                  !:
-       INTEGER(iwp)                               ::  m                  !:
-       INTEGER(iwp)                               ::  n                  !:
-       INTEGER(iwp)                               ::  npx                !:
-       INTEGER(iwp)                               ::  npy                !:
-       INTEGER(iwp)                               ::  nrx                !:
-       INTEGER(iwp)                               ::  nry                !:
-       INTEGER(iwp)                               ::  px                 !:
-       INTEGER(iwp)                               ::  py                 !:
-       INTEGER(iwp), DIMENSION(2)                 ::  scoord             !:
-       INTEGER(iwp)                               ::  server_pe          !:
-       INTEGER(iwp), DIMENSION(2)                 ::  size_of_array      !:
-
 
        IF ( myid == 0 )  THEN
+!--       TO_DO: Klaus: give more specific comment what size_of_array stands for
           CALL pmc_recv_from_client( client_id, size_of_array, 2, 0, 40, ierr )
           ALLOCATE( coarse_bound_all(size_of_array(1),size_of_array(2)) )
@@ -628,11 +648,19 @@
           CALL MPI_COMM_SIZE( comm1dx, npx, ierr )
           CALL MPI_COMM_SIZE( comm1dy, npy, ierr )
-
-          nrx = nxr - nxl + 1   !  +1 in index because FORTRAN indexing starts with 1, palm with 0
+!
+!--       The +1 in index is because PALM starts with nx=0
+          nrx = nxr - nxl + 1
           nry = nyn - nys + 1
           ic  = 0
-          DO  k = 1, size_of_array(2)                                   !  loop over all client PEs
-             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)       !  area in y required by actual client PE
-                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)    !  area in x required by actual client PE
+!
+!--       Loop over all client PEs
+          DO  k = 1, size_of_array(2)
+!
+!--          Area along y required by actual client PE
+             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
+!
+!--             Area along x required by actual client PE
+                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
+
                    px = i / nrx
                    py = j / nry
@@ -642,16 +670,34 @@
                   
                    ic = ic + 1
-                   index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp       !  First index in Server Array
-                   index_list(2,ic) = j - ( py * nry ) + 1 + nbgp       !  Second index in Server Array
-                   index_list(3,ic) = i - coarse_bound_all(1,k) + 1     !  x Index client coarse grid
-                   index_list(4,ic) = j - coarse_bound_all(3,k) + 1     !  y Index client coarse grid
-                   index_list(5,ic) = k - 1                             !  PE Number client
-                   index_list(6,ic) = server_pe                         !  PE Number server
+!
+!--                First index in server array
+                   index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp
+!
+!--                Second index in server array
+                   index_list(2,ic) = j - ( py * nry ) + 1 + nbgp
+!
+!--                x index of client coarse grid
+                   index_list(3,ic) = i - coarse_bound_all(1,k) + 1
+!
+!--                y index of client coarse grid
+                   index_list(4,ic) = j - coarse_bound_all(3,k) + 1
+!
+!--                PE number of client
+                   index_list(5,ic) = k - 1
+!
+!--                PE number of server
+                   index_list(6,ic) = server_pe
+
                 ENDDO
              ENDDO
           ENDDO
+!
+!--       TO_DO: Klaus: comment what is done here
           CALL pmc_s_set_2d_index_list( client_id, index_list(:,1:ic) )
+
        ELSE
-          ALLOCATE( index_list(6,1) )    !  Dummy allocate
+!
+!--       TO_DO: Klaus: comment why thie dummy allocation is required
+          ALLOCATE( index_list(6,1) )
           CALL pmc_s_set_2d_index_list( client_id, index_list )
        ENDIF
@@ -671,25 +717,28 @@
     IMPLICIT NONE
 
-    CHARACTER(LEN=DA_Namelen)  ::  myname     !:
-
-    INTEGER(iwp)               ::  i          !:
-    INTEGER(iwp)               ::  ierr       !:
-    INTEGER(iwp)               ::  icl        !:
-    INTEGER(iwp)               ::  icr        !:
-    INTEGER(iwp)               ::  j          !:
-    INTEGER(iwp)               ::  jcn        !:
-    INTEGER(iwp)               ::  jcs        !:
+    CHARACTER(LEN=da_namelen) ::  myname     !:
+
+    INTEGER(iwp) ::  i          !:
+    INTEGER(iwp) ::  ierr       !:
+    INTEGER(iwp) ::  icl        !:
+    INTEGER(iwp) ::  icr        !:
+    INTEGER(iwp) ::  j          !:
+    INTEGER(iwp) ::  jcn        !:
+    INTEGER(iwp) ::  jcs        !:
+
     INTEGER(iwp), DIMENSION(5) ::  val        !:
     
-    REAL(wp), DIMENSION(1)     ::  fval       !:
-    REAL(wp)                   ::  xcs        !:
-    REAL(wp)                   ::  xce        !:
-    REAL(wp)                   ::  ycs        !:
-    REAL(wp)                   ::  yce        !:
+    REAL(wp) ::  xcs        !:
+    REAL(wp) ::  xce        !:
+    REAL(wp) ::  ycs        !:
+    REAL(wp) ::  yce        !:
+
+    REAL(wp), DIMENSION(1) ::  fval       !:
                                              
-
+!
 !-- TO_DO: describe what is happening in this if-clause
 !-- Root Model does not have Server and is not a client
     IF ( .NOT. pmc_is_rootmodel() )  THEN
+
        CALL pmc_clientinit
 !
@@ -725,4 +774,5 @@
 
        IF ( myid == 0 )  THEN
+
           CALL pmc_send_to_server( val, SIZE( val ), 0, 123, ierr )
           CALL pmc_send_to_server( fval, SIZE( fval ), 0, 124, ierr )
@@ -732,4 +782,5 @@
 !
 !--       Receive Coarse grid information.
+!--       TO_DO: find shorter and more meaningful name for  define_coarse_grid_real
           CALL pmc_recv_from_server( define_coarse_grid_real,                  &
                                      SIZE(define_coarse_grid_real), 0, 21, ierr )
@@ -738,17 +789,18 @@
 !
 !--       Receive also the dz-,zu- and zw-arrays here.
-!--       TO_DO: what is the meaning of above comment + remove write statements
-!--              and give this informations in header
-          WRITE(0,*) 'Coarse grid from Server '
-          WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
-          WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
-          WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(5)
-          WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(6)
-          WRITE(0,*) 'dx            = ',define_coarse_grid_real(3)
-          WRITE(0,*) 'dy            = ',define_coarse_grid_real(4)
-          WRITE(0,*) 'dz            = ',define_coarse_grid_real(7)
-          WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
-          WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
-          WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)        
+!--       TO_DO: what is the meaning of above comment
+!
+!--        Debug-printouts - keep them
+!          WRITE(0,*) 'Coarse grid from Server '
+!          WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
+!          WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
+!          WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(5)
+!          WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(6)
+!          WRITE(0,*) 'dx            = ',define_coarse_grid_real(3)
+!          WRITE(0,*) 'dy            = ',define_coarse_grid_real(4)
+!          WRITE(0,*) 'dz            = ',define_coarse_grid_real(7)
+!          WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
+!          WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
+!          WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)
        ENDIF
 
@@ -765,5 +817,5 @@
 
 !
-!--    Get Server coordinates on coarse grid
+!--    Get server coordinates on coarse grid
        ALLOCATE( cg%coord_x(-nbgp:cg%nx+nbgp) )
        ALLOCATE( cg%coord_y(-nbgp:cg%ny+nbgp) )
@@ -776,31 +828,32 @@
 !
 !--    Get coarse grid coordinates and vales of the z-direction from server
-       IF ( myid == 0) THEN
-          CALL pmc_recv_from_server( cg%coord_x, cg%nx + 1 + 2 * nbgp, 0, 24,  &
-                                     ierr )
-          CALL pmc_recv_from_server( cg%coord_y, cg%ny + 1 + 2 * nbgp, 0, 25,  &
-                                     ierr )
+       IF ( myid == 0)  THEN
+
+          CALL pmc_recv_from_server( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr )
+          CALL pmc_recv_from_server( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr )
           CALL pmc_recv_from_server( cg%dzu, cg%nz + 1, 0, 26, ierr )
           CALL pmc_recv_from_server( cg%dzw, cg%nz + 1, 0, 27, ierr )
           CALL pmc_recv_from_server( cg%zu, cg%nz + 2, 0, 28, ierr )
           CALL pmc_recv_from_server( cg%zw, cg%nz + 2, 0, 29, ierr )
-       ENDIF
-
-!
-!--    and 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 )
+
+       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 )
        
+!
+!--    TO_DO: give comments what is happening here
        CALL pmci_map_fine_to_coarse_grid
        CALL pmc_c_get_2d_index_list
 
 !
-!--    Include couple arrays into client content.
+!--    Include couple arrays into client content
+!--    TO_DO: Klaus: better explain the above comment (what is client content?)
        CALL  pmc_c_clear_next_array_list
        DO  WHILE ( pmc_c_getnextarray( myname ) )
@@ -823,5 +876,6 @@
 
 !
-!--    Two-way coupling      
+!--    Two-way coupling
+!--    TO_DO: comment what is happening here
        IF ( nesting_mode == 'two-way' )  THEN
           CALL pmci_init_anterp_tophat
@@ -837,507 +891,507 @@
  CONTAINS
 
-
     SUBROUTINE pmci_map_fine_to_coarse_grid
 
-        IMPLICIT NONE
-
-        INTEGER(iwp), DIMENSION(5,numprocs)  ::  coarse_bound_all   !:
-        INTEGER(iwp), DIMENSION(2)           ::  size_of_array      !:
+       IMPLICIT NONE
+
+       INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all   !:
+       INTEGER(iwp), DIMENSION(2)          ::  size_of_array      !:
                                               
-        REAL(wp)                             ::  coarse_dx   !:
-        REAL(wp)                             ::  coarse_dy   !:
-        REAL(wp)                             ::  loffset     !:
-        REAL(wp)                             ::  noffset     !:
-        REAL(wp)                             ::  roffset     !:
-        REAL(wp)                             ::  soffset     !:
-
-!
-!--     Determine indices of interpolation/anterpolation area in coarse grid.
-        coarse_dx = cg%dx
-        coarse_dy = cg%dy
-        
-        loffset = MOD( coord_x(nxl), coarse_dx )       !  If the fine- and coarse grid nodes do not match.
-        xexl    = coarse_dx + loffset
-        nhll    = CEILING( xexl / coarse_dx )          !  This is needed in the anterpolation phase.
-        xcs     = coord_x(nxl) - xexl
-        DO  i = 0, cg%nx
-           IF ( cg%coord_x(i) > xcs )  THEN
-              icl = MAX( -1, i-1 )  
-              EXIT
-           ENDIF
-        ENDDO
-
-        roffset = MOD( coord_x(nxr + 1), coarse_dx )   !  If the fine- and coarse grid nodes do not match.
-        xexr    = coarse_dx + roffset
-        nhlr    = CEILING( xexr / coarse_dx )          !  This is needed in the anterpolation phase.
-        xce     = coord_x(nxr) + xexr
-        DO  i = cg%nx, 0 , -1
-           IF ( cg%coord_x(i) < xce )  THEN
-              icr = MIN( cg%nx + 1, i + 1 )
-              EXIT
-           ENDIF
-        ENDDO
-
-        soffset = MOD( coord_y(nys), coarse_dy )       !  If the fine- and coarse grid nodes do not match
-        yexs    = coarse_dy + soffset
-        nhls    = CEILING( yexs / coarse_dy )          !  This is needed in the anterpolation phase.
-        ycs     = coord_y(nys) - yexs
-        DO  j = 0, cg%ny
-           IF ( cg%coord_y(j) > ycs )  THEN
-              jcs = MAX( -nbgp, j - 1 )  
-              EXIT
-           ENDIF
-        ENDDO
-
-        noffset = MOD( coord_y(nyn + 1), coarse_dy )   !  If the fine- and coarse grid nodes do not match
-        yexn    = coarse_dy + noffset
-        nhln    = CEILING( yexn / coarse_dy )          !  This is needed in the anterpolation phase.
-        yce     = coord_y(nyn) + yexn
-        DO  j = cg%ny, 0, -1
-           IF ( cg%coord_y(j) < yce )  THEN
-              jcn = MIN( cg%ny + nbgp, j + 1 )  
-              EXIT
-           ENDIF
-        ENDDO
-
-        coarse_bound(1) = icl
-        coarse_bound(2) = icr
-        coarse_bound(3) = jcs
-        coarse_bound(4) = jcn
-        coarse_bound(5) = myid
-!
-!--     Note that MPI_Gather receives data from all processes in the rank order. 
-        CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5, &
-                         MPI_INTEGER, 0, comm2d, ierr )
-
-        IF ( myid == 0 )  THEN
-           size_of_array(1) = SIZE( coarse_bound_all, 1 )
-           size_of_array(2) = SIZE( coarse_bound_all, 2 )
-           CALL pmc_send_to_server( size_of_array, 2, 0, 40, ierr )
-           CALL pmc_send_to_server( coarse_bound_all, SIZE( coarse_bound_all ), 0, 41, ierr )
-        ENDIF
-
-   END SUBROUTINE pmci_map_fine_to_coarse_grid
-
-
-
-   SUBROUTINE pmci_init_interp_tril
-
-!
-!--   Precomputation of the interpolation coefficients and client-array indices 
-!--   to be used by the interpolation routines interp_tril_lr, interp_tril_ns and 
-!--   interp_tril_t. Constant dz is still assumed.
-!
-!--                                  Antti Hellsten 3.3.2015.
-!
-!--   Modified for variable dz, but not yet tested.
-!--                                  Antti Hellsten 27.3.2015.
-!
-      IMPLICIT NONE
-
-      INTEGER(iwp) ::  i       !:
-      INTEGER(iwp) ::  i1      !:
-      INTEGER(iwp) ::  j       !:
-      INTEGER(iwp) ::  j1      !:
-      INTEGER(iwp) ::  k       !:
-      INTEGER(iwp) ::  kc      !:
-
-      REAL(wp) ::  coarse_dx   !:
-      REAL(wp) ::  coarse_dy   !:
-      REAL(wp) ::  coarse_dz   !:
-      REAL(wp) ::  xb          !:
-      REAL(wp) ::  xcsu        !:
-      REAL(wp) ::  xfso        !:
-      REAL(wp) ::  xcso        !:
-      REAL(wp) ::  xfsu        !:
-      REAL(wp) ::  yb          !:
-      REAL(wp) ::  ycso        !:
-      REAL(wp) ::  ycsv        !:
-      REAL(wp) ::  yfso        !:
-      REAL(wp) ::  yfsv        !:
-      REAL(wp) ::  zcso        !:
-      REAL(wp) ::  zcsw        !:
-      REAL(wp) ::  zfso        !:
-      REAL(wp) ::  zfsw        !:
+       REAL(wp) ::  loffset     !:
+       REAL(wp) ::  noffset     !:
+       REAL(wp) ::  roffset     !:
+       REAL(wp) ::  soffset     !:
+
+!
+!--    Determine indices of interpolation/anterpolation area in the coarse grid
+!--    If the fine- and coarse grid nodes do not match.
+       loffset = MOD( coord_x(nxl), cg%dx )
+       xexl    = cg%dx + loffset
+!
+!--    This is needed in the anterpolation phase
+       nhll = CEILING( xexl / cg%dx )
+       xcs  = coord_x(nxl) - xexl
+       DO  i = 0, cg%nx
+          IF ( cg%coord_x(i) > xcs )  THEN
+             icl = MAX( -1, i-1 )
+             EXIT
+          ENDIF
+       ENDDO
+!
+!--    If the fine- and coarse grid nodes do not match
+       roffset = MOD( coord_x(nxr+1), cg%dx )
+       xexr    = cg%dx + roffset
+!
+!--    This is needed in the anterpolation phase
+       nhlr = CEILING( xexr / cg%dx )
+       xce  = coord_x(nxr) + xexr
+       DO  i = cg%nx, 0 , -1
+          IF ( cg%coord_x(i) < xce )  THEN
+             icr = MIN( cg%nx+1, i+1 )
+             EXIT
+          ENDIF
+       ENDDO
+!
+!--    If the fine- and coarse grid nodes do not match
+       soffset = MOD( coord_y(nys), cg%dy )
+       yexs    = cg%dy + soffset
+!
+!--    This is needed in the anterpolation phase
+       nhls = CEILING( yexs / cg%dy )
+       ycs  = coord_y(nys) - yexs
+       DO  j = 0, cg%ny
+          IF ( cg%coord_y(j) > ycs )  THEN
+             jcs = MAX( -nbgp, j-1 )
+             EXIT
+          ENDIF
+       ENDDO
+!
+!--    If the fine- and coarse grid nodes do not match
+       noffset = MOD( coord_y(nyn+1), cg%dy )
+       yexn    = cg%dy + noffset
+!
+!--    This is needed in the anterpolation phase
+       nhln = CEILING( yexn / cg%dy )
+       yce  = coord_y(nyn) + yexn
+       DO  j = cg%ny, 0, -1
+          IF ( cg%coord_y(j) < yce )  THEN
+             jcn = MIN( cg%ny + nbgp, j+1 )
+             EXIT
+          ENDIF
+       ENDDO
+
+       coarse_bound(1) = icl
+       coarse_bound(2) = icr
+       coarse_bound(3) = jcs
+       coarse_bound(4) = jcn
+       coarse_bound(5) = myid
+!
+!--    Note that MPI_Gather receives data from all processes in the rank order
+!--    TO_DO: refer to the line where this fact becomes important
+       CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5, &
+                        MPI_INTEGER, 0, comm2d, ierr )
+
+       IF ( myid == 0 )  THEN
+          size_of_array(1) = SIZE( coarse_bound_all, 1 )
+          size_of_array(2) = SIZE( coarse_bound_all, 2 )
+          CALL pmc_send_to_server( size_of_array, 2, 0, 40, ierr )
+          CALL pmc_send_to_server( coarse_bound_all, SIZE( coarse_bound_all ), &
+                                   0, 41, ierr )
+       ENDIF
+
+    END SUBROUTINE pmci_map_fine_to_coarse_grid
+
+
+
+    SUBROUTINE pmci_init_interp_tril
+!
+!--    Precomputation of the interpolation coefficients and client-array indices
+!--    to be used by the interpolation routines interp_tril_lr, interp_tril_ns
+!--    and interp_tril_t. Constant dz is still assumed.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i       !:
+       INTEGER(iwp) ::  i1      !:
+       INTEGER(iwp) ::  j       !:
+       INTEGER(iwp) ::  j1      !:
+       INTEGER(iwp) ::  k       !:
+       INTEGER(iwp) ::  kc      !:
+
+       REAL(wp) ::  xb          !:
+       REAL(wp) ::  xcsu        !:
+       REAL(wp) ::  xfso        !:
+       REAL(wp) ::  xcso        !:
+       REAL(wp) ::  xfsu        !:
+       REAL(wp) ::  yb          !:
+       REAL(wp) ::  ycso        !:
+       REAL(wp) ::  ycsv        !:
+       REAL(wp) ::  yfso        !:
+       REAL(wp) ::  yfsv        !:
+       REAL(wp) ::  zcso        !:
+       REAL(wp) ::  zcsw        !:
+       REAL(wp) ::  zfso        !:
+       REAL(wp) ::  zfsw        !:
       
 
-      coarse_dx = cg%dx
-      coarse_dy = cg%dy
-      coarse_dz = cg%dz
-      xb        = nxl * dx 
-      yb        = nys * dy
+       xb = nxl * dx
+       yb = nys * dy
       
-      ALLOCATE( icu(nxlg:nxrg) )
-      ALLOCATE( ico(nxlg:nxrg) )
-      ALLOCATE( jcv(nysg:nyng) )
-      ALLOCATE( jco(nysg:nyng) )
-      ALLOCATE( kcw(nzb:nzt+1) )
-      ALLOCATE( kco(nzb:nzt+1) )
-      ALLOCATE( r1xu(nxlg:nxrg) )
-      ALLOCATE( r2xu(nxlg:nxrg) )
-      ALLOCATE( r1xo(nxlg:nxrg) )
-      ALLOCATE( r2xo(nxlg:nxrg) )
-      ALLOCATE( r1yv(nysg:nyng) )
-      ALLOCATE( r2yv(nysg:nyng) )
-      ALLOCATE( r1yo(nysg:nyng) )
-      ALLOCATE( r2yo(nysg:nyng) )
-      ALLOCATE( r1zw(nzb:nzt+1) )
-      ALLOCATE( r2zw(nzb:nzt+1) )
-      ALLOCATE( r1zo(nzb:nzt+1) )
-      ALLOCATE( r2zo(nzb:nzt+1) )
-
-!
-!--   Note that the node coordinates xfs... and xcs... are relative 
-!--   to the lower-left-bottom corner of the fc-array, not the actual 
-!--   client domain corner.
-      DO  i = nxlg, nxrg
-         xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
-         xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
-         icu(i)  = icl + FLOOR( xfsu / coarse_dx )
-         ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * coarse_dx ) / coarse_dx )
-         xcsu    = ( icu(i) - icl ) * coarse_dx 
-         xcso    = ( ico(i) - icl ) * coarse_dx + 0.5_wp * coarse_dx
-         r2xu(i) = ( xfsu - xcsu ) / coarse_dx
-         r2xo(i) = ( xfso - xcso ) / coarse_dx
-         r1xu(i) = 1.0_wp - r2xu(i)
-         r1xo(i) = 1.0_wp - r2xo(i)
-      ENDDO
-
-      DO  j = nysg, nyng
-         yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
-         yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
-         jcv(j)  = jcs + FLOOR( yfsv/coarse_dy )           
-         jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * coarse_dy ) / coarse_dy )
-         ycsv    = ( jcv(j) - jcs ) * coarse_dy
-         ycso    = ( jco(j) - jcs ) * coarse_dy + 0.5_wp * coarse_dy         
-         r2yv(j) = ( yfsv - ycsv ) / coarse_dy            
-         r2yo(j) = ( yfso - ycso ) / coarse_dy  
-         r1yv(j) = 1.0_wp - r2yv(j)
-         r1yo(j) = 1.0_wp - r2yo(j)        
-      ENDDO
-
-      DO  k = nzb, nzt + 1
-         zfsw = zw(k)
-         zfso = zu(k)
-
-         kc = 0
-         DO  WHILE ( cg%zw(kc) <= zfsw ) 
-            kc = kc + 1
-         ENDDO
-         kcw(k) = kc - 1
+       ALLOCATE( icu(nxlg:nxrg) )
+       ALLOCATE( ico(nxlg:nxrg) )
+       ALLOCATE( jcv(nysg:nyng) )
+       ALLOCATE( jco(nysg:nyng) )
+       ALLOCATE( kcw(nzb:nzt+1) )
+       ALLOCATE( kco(nzb:nzt+1) )
+       ALLOCATE( r1xu(nxlg:nxrg) )
+       ALLOCATE( r2xu(nxlg:nxrg) )
+       ALLOCATE( r1xo(nxlg:nxrg) )
+       ALLOCATE( r2xo(nxlg:nxrg) )
+       ALLOCATE( r1yv(nysg:nyng) )
+       ALLOCATE( r2yv(nysg:nyng) )
+       ALLOCATE( r1yo(nysg:nyng) )
+       ALLOCATE( r2yo(nysg:nyng) )
+       ALLOCATE( r1zw(nzb:nzt+1) )
+       ALLOCATE( r2zw(nzb:nzt+1) )
+       ALLOCATE( r1zo(nzb:nzt+1) )
+       ALLOCATE( r2zo(nzb:nzt+1) )
+
+!
+!--    Note that the node coordinates xfs... and xcs... are relative to the
+!--    lower-left-bottom corner of the fc-array, not the actual client domain
+!--    corner
+       DO  i = nxlg, nxrg
+          xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
+          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
+          icu(i)  = icl + FLOOR( xfsu / cg%dx )
+          ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * cg%dx ) / cg%dx )
+          xcsu    = ( icu(i) - icl ) * cg%dx
+          xcso    = ( ico(i) - icl ) * cg%dx + 0.5_wp * cg%dx
+          r2xu(i) = ( xfsu - xcsu ) / cg%dx
+          r2xo(i) = ( xfso - xcso ) / cg%dx
+          r1xu(i) = 1.0_wp - r2xu(i)
+          r1xo(i) = 1.0_wp - r2xo(i)
+       ENDDO
+
+       DO  j = nysg, nyng
+          yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
+          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
+          jcv(j)  = jcs + FLOOR( yfsv / cg%dy )
+          jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * cg%dy ) / cg%dy )
+          ycsv    = ( jcv(j) - jcs ) * cg%dy
+          ycso    = ( jco(j) - jcs ) * cg%dy + 0.5_wp * cg%dy
+          r2yv(j) = ( yfsv - ycsv ) / cg%dy
+          r2yo(j) = ( yfso - ycso ) / cg%dy
+          r1yv(j) = 1.0_wp - r2yv(j)
+          r1yo(j) = 1.0_wp - r2yo(j)
+       ENDDO
+
+       DO  k = nzb, nzt + 1
+          zfsw = zw(k)
+          zfso = zu(k)
+
+          kc = 0
+          DO  WHILE ( cg%zw(kc) <= zfsw )
+             kc = kc + 1
+          ENDDO
+          kcw(k) = kc - 1
          
-         kc = 0
-         DO  WHILE ( cg%zu(kc) <= zfso ) 
-            kc = kc + 1
-         ENDDO
-         kco(k) = kc - 1
-
-         zcsw    = cg%zw(kcw(k))
-         zcso    = cg%zu(kco(k))
-         r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k) + 1 )
-         r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k) + 1 )
-         r1zw(k) = 1.0_wp - r2zw(k)
-         r1zo(k) = 1.0_wp - r2zo(k)
-      ENDDO
+          kc = 0
+          DO  WHILE ( cg%zu(kc) <= zfso )
+             kc = kc + 1
+          ENDDO
+          kco(k) = kc - 1
+
+          zcsw    = cg%zw(kcw(k))
+          zcso    = cg%zu(kco(k))
+          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k)+1)
+          r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k)+1)
+          r1zw(k) = 1.0_wp - r2zw(k)
+          r1zo(k) = 1.0_wp - r2zo(k)
+       ENDDO
       
-   END SUBROUTINE pmci_init_interp_tril
-
-
-
-   SUBROUTINE pmci_init_loglaw_correction
-
-!
-!--   Precomputation of the index and log-ratio arrays for the log-law corrections
-!--   for near-wall nodes after the nest-BC interpolation.
-!--   These are used by the interpolation routines interp_tril_lr and interp_tril_ns.
-!
-!--                                  Antti Hellsten 30.12.2015.
-!
-      IMPLICIT NONE
-      INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i.  
-      INTEGER(iwp) ::  i            !:
-      INTEGER(iwp) ::  icorr        !:
-      INTEGER(iwp) ::  inc          !:  Wall outward-normal index increment -1 or 1, for direction=1, inc=1 always. 
-      INTEGER(iwp) ::  iw           !:
-      INTEGER(iwp) ::  j            !:
-      INTEGER(iwp) ::  jcorr        !:
-      INTEGER(iwp) ::  jw           !:
-      INTEGER(iwp) ::  k            !:
-      INTEGER(iwp) ::  kb           !:
-      INTEGER(iwp) ::  kcorr        !:
-      INTEGER(iwp) ::  lc           !:
-      INTEGER(iwp) ::  ni           !:
-      INTEGER(iwp) ::  nj           !:
-      INTEGER(iwp) ::  nk           !:
-      INTEGER(iwp) ::  nzt_topo_max !:
-      INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
-
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
-
-!
-!--   First determine the maximum k-index needed for the near-wall corrections.
-!--   This maximum is individual for each boundary to minimize the storage requirements
-!--   and to minimize the corresponding loop k-range in the interpolation routines.
-      nzt_topo_nestbc_l = nzb
-      IF ( nest_bound_l )  THEN
-         DO  i = nxl - 1, nxl
-            DO  j = nys, nyn
-               nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),   &
-                                        nzb_v_inner(j,i), nzb_w_inner(j,i) )
-            ENDDO
-         ENDDO
-         nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
-      ENDIF
+    END SUBROUTINE pmci_init_interp_tril
+
+
+
+    SUBROUTINE pmci_init_loglaw_correction
+!
+!--    Precomputation of the index and log-ratio arrays for the log-law
+!--    corrections for near-wall nodes after the nest-BC interpolation.
+!--    These are used by the interpolation routines interp_tril_lr and
+!--    interp_tril_ns.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i.
+       INTEGER(iwp) ::  i            !:
+       INTEGER(iwp) ::  icorr        !:
+       INTEGER(iwp) ::  inc          !:  Wall outward-normal index increment -1
+                                     !: or 1, for direction=1, inc=1 always
+       INTEGER(iwp) ::  iw           !:
+       INTEGER(iwp) ::  j            !:
+       INTEGER(iwp) ::  jcorr        !:
+       INTEGER(iwp) ::  jw           !:
+       INTEGER(iwp) ::  k            !:
+       INTEGER(iwp) ::  kb           !:
+       INTEGER(iwp) ::  kcorr        !:
+       INTEGER(iwp) ::  lc           !:
+       INTEGER(iwp) ::  ni           !:
+       INTEGER(iwp) ::  nj           !:
+       INTEGER(iwp) ::  nk           !:
+       INTEGER(iwp) ::  nzt_topo_max !:
+       INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
+
+       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
+
+!
+!--    First determine the maximum k-index needed for the near-wall corrections.
+!--    This maximum is individual for each boundary to minimize the storage
+!--    requirements and to minimize the corresponding loop k-range in the
+!--    interpolation routines.
+       nzt_topo_nestbc_l = nzb
+       IF ( nest_bound_l )  THEN
+          DO  i = nxl-1, nxl
+             DO  j = nys, nyn
+                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),  &
+                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
+             ENDDO
+          ENDDO
+          nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
+       ENDIF
       
-      nzt_topo_nestbc_r = nzb
-      IF ( nest_bound_r )  THEN
-         i = nxr + 1
-         DO  j = nys, nyn
-            nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),      &
-                                     nzb_v_inner(j,i), nzb_w_inner(j,i) )
-         ENDDO
-         nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1       
-      ENDIF
-
-      nzt_topo_nestbc_s = nzb
-      IF ( nest_bound_s )  THEN
-         DO  j = nys - 1, nys
-            DO  i = nxl, nxr
-               nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),   &
-                                        nzb_v_inner(j,i), nzb_w_inner(j,i) )
-            ENDDO
-         ENDDO
-         nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1      
-      ENDIF
-
-      nzt_topo_nestbc_n = nzb
-      IF ( nest_bound_n )  THEN
-         j = nyn + 1
-         DO  i = nxl, nxr
-            nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),      &
-                                     nzb_v_inner(j,i), nzb_w_inner(j,i) )
-         ENDDO
-         nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
-      ENDIF
-
-!
-!--   Then determine the maximum number of near-wall nodes per wall point based
-!--   on the grid-spacing ratios.
-      nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,                &
-                          nzt_topo_nestbc_s, nzt_topo_nestbc_n )
-
-!
-!--   Note that the outer division must be integer division.
-      ni = CEILING( cg%dx / dx ) / 2
-      nj = CEILING( cg%dy / dy ) / 2
-      nk = 1
-      DO  k = 1, nzt_topo_max
-         nk = MAX( nk, CEILING( cg%dzu(k) / dzu(k) ) )
-      ENDDO
-      nk = nk / 2   !  Note that this must be integer division.
-      ncorr =  MAX( ni, nj, nk )
-
-      ALLOCATE( lcr(0:ncorr - 1) )
-      lcr = 1.0_wp
-
-!
-!--   First horizontal walls
-!--   Left boundary
-      IF ( nest_bound_l )   THEN
-         ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
-         ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
-         ALLOCATE( logc_ratio_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
-         ALLOCATE( logc_ratio_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
-         logc_u_l       = 0
-         logc_v_l       = 0
-         logc_ratio_u_l = 1.0_wp
-         logc_ratio_v_l = 1.0_wp
-         direction      = 1
-         inc            = 1
-
-         DO  j = nys, nyn      
-!
-!--         Left boundary for u
-            i   = 0
-            kb  = nzb_u_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_u_l(k,j,1) = lc
-            logc_ratio_u_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-!
-!--         Left boundary for v
-            i   = -1
-            kb  =  nzb_v_inner(j,i)
-            k   =  kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_v_l(k,j,1) = lc
-            logc_ratio_v_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-         ENDDO
-      ENDIF
-
-!
-!--   Right boundary
-      IF ( nest_bound_r )  THEN
-         ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
-         ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
-         ALLOCATE( logc_ratio_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
-         ALLOCATE( logc_ratio_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
-         logc_u_r       = 0
-         logc_v_r       = 0
-         logc_ratio_u_r = 1.0_wp
-         logc_ratio_v_r = 1.0_wp
-         direction      = 1
-         inc            = 1  
-         DO  j = nys, nyn
-!
-!--         Right boundary for u.
-            i   = nxr + 1
-            kb  = nzb_u_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_u_r(k,j,1) = lc
-            logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-!
-!--         Right boundary for v.
-            i   = nxr + 1
-            kb  = nzb_v_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_v_r(k,j,1) = lc
-            logc_ratio_v_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-         ENDDO
-      ENDIF
-
-!
-!--   South boundary
-      IF ( nest_bound_s )  THEN
-         ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
-         ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
-         ALLOCATE( logc_ratio_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
-         ALLOCATE( logc_ratio_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
-         logc_u_s       = 0
-         logc_v_s       = 0
-         logc_ratio_u_s = 1.0_wp
-         logc_ratio_v_s = 1.0_wp
-         direction      = 1
-         inc            = 1 
-         DO  i = nxl, nxr
-!
-!--         South boundary for u.
-            j   = -1
-            kb  =  nzb_u_inner(j,i)
-            k   =  kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_u_s(k,i,1) = lc
-            logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-!
-!--         South boundary for v
-            j   = 0
-            kb  = nzb_v_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_v_s(k,i,1) = lc
-            logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-         ENDDO
-      ENDIF
-
-!
-!--   North boundary
-      IF ( nest_bound_n )  THEN
-         ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
-         ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
-         ALLOCATE( logc_ratio_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
-         ALLOCATE( logc_ratio_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
-         logc_u_n       = 0
-         logc_v_n       = 0
-         logc_ratio_u_n = 1.0_wp
-         logc_ratio_v_n = 1.0_wp 
-         direction      = 1
-         inc            = 1 
-         DO  i = nxl, nxr
-!
-!--         North boundary for u.
-            j   = nyn + 1
-            kb  = nzb_u_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_u_n(k,i,1) = lc
-            logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-!
-!--         North boundary for v.
-            j   = nyn + 1
-            kb  = nzb_v_inner(j,i)
-            k   = kb + 1
-            wall_index = kb 
-            CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
-                                     wall_index, z0(j,i), kb, direction, ncorr )
-            logc_v_n(k,i,1) = lc
-            logc_ratio_v_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
-            lcr(0:ncorr-1) = 1.0_wp
-
-         ENDDO
-      ENDIF
+       nzt_topo_nestbc_r = nzb
+       IF ( nest_bound_r )  THEN
+          i = nxr + 1
+          DO  j = nys, nyn
+             nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),     &
+                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
+          ENDDO
+          nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
+       ENDIF
+
+       nzt_topo_nestbc_s = nzb
+       IF ( nest_bound_s )  THEN
+          DO  j = nys-1, nys
+             DO  i = nxl, nxr
+                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),  &
+                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
+             ENDDO
+          ENDDO
+          nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1
+       ENDIF
+
+       nzt_topo_nestbc_n = nzb
+       IF ( nest_bound_n )  THEN
+          j = nyn + 1
+          DO  i = nxl, nxr
+             nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),     &
+                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
+          ENDDO
+          nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
+       ENDIF
+
+!
+!--    Then determine the maximum number of near-wall nodes per wall point based
+!--    on the grid-spacing ratios.
+       nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,               &
+                           nzt_topo_nestbc_s, nzt_topo_nestbc_n )
+
+!
+!--    Note that the outer division must be integer division.
+       ni = CEILING( cg%dx / dx ) / 2
+       nj = CEILING( cg%dy / dy ) / 2
+       nk = 1
+       DO  k = 1, nzt_topo_max
+          nk = MAX( nk, CEILING( cg%dzu(k) / dzu(k) ) )
+       ENDDO
+       nk = nk / 2   !  Note that this must be integer division.
+       ncorr =  MAX( ni, nj, nk )
+
+       ALLOCATE( lcr(0:ncorr-1) )
+       lcr = 1.0_wp
+
+!
+!--    First horizontal walls
+!--    Left boundary
+       IF ( nest_bound_l )  THEN
+          ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
+          ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
+          ALLOCATE( logc_ratio_u_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
+          ALLOCATE( logc_ratio_v_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
+          logc_u_l       = 0
+          logc_v_l       = 0
+          logc_ratio_u_l = 1.0_wp
+          logc_ratio_v_l = 1.0_wp
+          direction      = 1
+          inc            = 1
+
+          DO  j = nys, nyn
+!
+!--          Left boundary for u
+             i   = 0
+             kb  = nzb_u_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_u_l(k,j,1) = lc
+             logc_ratio_u_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+!
+!--          Left boundary for v
+             i   = -1
+             kb  =  nzb_v_inner(j,i)
+             k   =  kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_v_l(k,j,1) = lc
+             logc_ratio_v_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+          ENDDO
+       ENDIF
+
+!
+!--    Right boundary
+       IF ( nest_bound_r )  THEN
+          ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
+          ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
+          ALLOCATE( logc_ratio_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
+          ALLOCATE( logc_ratio_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
+          logc_u_r       = 0
+          logc_v_r       = 0
+          logc_ratio_u_r = 1.0_wp
+          logc_ratio_v_r = 1.0_wp
+          direction      = 1
+          inc            = 1
+          DO  j = nys, nyn
+!
+!--          Right boundary for u.
+             i   = nxr + 1
+             kb  = nzb_u_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_u_r(k,j,1) = lc
+             logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+!
+!--          Right boundary for v.
+             i   = nxr + 1
+             kb  = nzb_v_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_v_r(k,j,1) = lc
+             logc_ratio_v_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+          ENDDO
+       ENDIF
+
+!
+!--    South boundary
+       IF ( nest_bound_s )  THEN
+          ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
+          ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
+          ALLOCATE( logc_ratio_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
+          ALLOCATE( logc_ratio_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
+          logc_u_s       = 0
+          logc_v_s       = 0
+          logc_ratio_u_s = 1.0_wp
+          logc_ratio_v_s = 1.0_wp
+          direction      = 1
+          inc            = 1
+          DO  i = nxl, nxr
+!
+!--          South boundary for u.
+             j   = -1
+             kb  =  nzb_u_inner(j,i)
+             k   =  kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_u_s(k,i,1) = lc
+             logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+!
+!--          South boundary for v
+             j   = 0
+             kb  = nzb_v_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_v_s(k,i,1) = lc
+             logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+          ENDDO
+       ENDIF
+
+!
+!--    North boundary
+       IF ( nest_bound_n )  THEN
+          ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
+          ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
+          ALLOCATE( logc_ratio_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
+          ALLOCATE( logc_ratio_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
+          logc_u_n       = 0
+          logc_v_n       = 0
+          logc_ratio_u_n = 1.0_wp
+          logc_ratio_v_n = 1.0_wp
+          direction      = 1
+          inc            = 1
+          DO  i = nxl, nxr
+!
+!--          North boundary for u.
+             j   = nyn + 1
+             kb  = nzb_u_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_u_n(k,i,1) = lc
+             logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+!
+!--          North boundary for v.
+             j   = nyn + 1
+             kb  = nzb_v_inner(j,i)
+             k   = kb + 1
+             wall_index = kb
+             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
+                                inc, wall_index, z0(j,i), kb, direction, ncorr )
+             logc_v_n(k,i,1) = lc
+             logc_ratio_v_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
+             lcr(0:ncorr-1) = 1.0_wp
+
+          ENDDO
+       ENDIF
 
 !       
-!--   Then vertical walls and corners if necessary.
-      IF ( topography /= 'flat' ) THEN        
-         kb = 0       ! kb is not used when direction > 1
+!--    Then vertical walls and corners if necessary.
+       IF ( topography /= 'flat' )  THEN
+          kb = 0       ! kb is not used when direction > 1
 !        
-!--      Left boundary
-         IF ( nest_bound_l )  THEN
-            ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
-            ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
-            logc_w_l       = 0                   
-            logc_ratio_w_l = 1.0_wp
-            direction      = 2
-            DO  j = nys, nyn
-               DO  k = nzb, nzt_topo_nestbc_l
-
-!
-!--               Wall for u on the south side, but not on the north side.
-                  i  = 0
-                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.     &
-                       ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )  THEN
-                     inc        =  1
-                     wall_index =  j 
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
-
-!
-!--                  The direction of the wall-normal index is stored as the sign of the logc-element. 
-                     logc_u_l(k,j,2) = inc * lc
-                     logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for u on the north side, but not on the south side.
-                  i  = 0
+!--       Left boundary
+          IF ( nest_bound_l )  THEN
+             ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
+             ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,       &
+                                      0:ncorr-1) )
+             logc_w_l       = 0
+             logc_ratio_w_l = 1.0_wp
+             direction      = 2
+             DO  j = nys, nyn
+                DO  k = nzb, nzt_topo_nestbc_l
+
+!
+!--                Wall for u on the south side, but not on the north side
+                   i  = 0
+                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.        &
+                        ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )           &
+                   THEN
+                      inc        =  1
+                      wall_index =  j
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
+!
+!--                   The direction of the wall-normal index is stored as the
+!--                   sign of the logc-element.
+                      logc_u_l(k,j,2) = inc * lc
+                      logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+!
+!--                Wall for u on the north side, but not on the south side
+                   i  = 0
+!--                TO_DO: routine must be indentet by 1 space from here on,
+!--                       and long lines must be wrapped
                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND.     &
                        ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
@@ -1619,211 +1673,230 @@
 
 
-   SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc, wall_index, z0_l, kb, direction, ncorr )
-      IMPLICIT NONE
-      INTEGER(iwp), INTENT(IN)  ::  direction                 !:
-      INTEGER(iwp), INTENT(IN)  ::  ij                        !:
-      INTEGER(iwp), INTENT(IN)  ::  inc                       !:
-      INTEGER(iwp), INTENT(IN)  ::  k                         !:
-      INTEGER(iwp), INTENT(IN)  ::  kb                        !:
-      INTEGER(iwp), INTENT(OUT) ::  lc                        !:
-      INTEGER(iwp), INTENT(IN)  ::  ncorr                     !:
-      INTEGER(iwp), INTENT(IN)  ::  wall_index                !:
-      REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !:
-      REAL(wp), INTENT(IN)      ::  z0_l                      !:
+    SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc,  &
+                                        wall_index, z0_l, kb, direction, ncorr )
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN)  ::  direction                 !:
+       INTEGER(iwp), INTENT(IN)  ::  ij                        !:
+       INTEGER(iwp), INTENT(IN)  ::  inc                       !:
+       INTEGER(iwp), INTENT(IN)  ::  k                         !:
+       INTEGER(iwp), INTENT(IN)  ::  kb                        !:
+       INTEGER(iwp), INTENT(OUT) ::  lc                        !:
+       INTEGER(iwp), INTENT(IN)  ::  ncorr                     !:
+       INTEGER(iwp), INTENT(IN)  ::  wall_index                !:
+
+       INTEGER(iwp) ::  alcorr       !:
+       INTEGER(iwp) ::  corr_index   !:
+       INTEGER(iwp) ::  lcorr        !:
+
+       LOGICAL      ::  more         !:
+
+       REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !:
+       REAL(wp), INTENT(IN)      ::  z0_l                      !:
       
-      INTEGER(iwp) ::  alcorr       !:
-      INTEGER(iwp) ::  corr_index   !:
-      INTEGER(iwp) ::  lcorr        !:
-      REAL(wp)     ::  logvelc1     !:
-      LOGICAL      ::  more         !:
+       REAL(wp)     ::  logvelc1     !:
       
 
-      SELECT CASE ( direction )
-
-         CASE (1)   !  k
-            more = .TRUE.
-            lcorr = 0
-            DO  WHILE ( more .AND. lcorr <= ncorr - 1 )
-               corr_index = k + lcorr 
-               IF ( lcorr == 0 )  THEN
-                  CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
-               ENDIF
+       SELECT CASE ( direction )
+
+          CASE (1)   !  k
+             more = .TRUE.
+             lcorr = 0
+             DO  WHILE ( more .AND. lcorr <= ncorr-1 )
+                corr_index = k + lcorr
+                IF ( lcorr == 0 )  THEN
+                   CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
+                ENDIF
                
-               IF ( corr_index < lc )  THEN
-                  lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
-                  more = .TRUE.
-               ELSE
-                  lcr(lcorr) = 1.0
-                  more = .FALSE.
-               ENDIF
-               lcorr = lcorr + 1              
-            ENDDO
-
-         CASE (2)   !  j 
-            more = .TRUE.
-            lcorr  = 0
-            alcorr = 0
-            DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
-               corr_index = ij + lcorr   !  In this case (direction = 2) ij is j
-               IF ( lcorr == 0 )  THEN
-                  CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index, z0_l, inc )
-               ENDIF
-
-!
-!--            The role of inc here is to make the comparison operation "<" valid in both directions.
-               IF ( inc * corr_index < inc * lc )  THEN 
-                  lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy  &
-                              - coord_y(wall_index) ) / z0_l ) / logvelc1
-                  more = .TRUE.
-               ELSE
-                  lcr(alcorr) = 1.0_wp
-                  more = .FALSE.
-               ENDIF
-               lcorr  = lcorr + inc
-               alcorr = ABS( lcorr )
-            ENDDO
-
-         CASE (3)   !  i          
-            more = .TRUE.
-            lcorr  = 0
-            alcorr = 0
-            DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
-               corr_index = ij + lcorr   !  In this case (direction = 3) ij is i
-               IF ( lcorr == 0 )  THEN
-                  CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index, z0_l, inc )
-               ENDIF
-
-!
-!--            The role of inc here is to make the comparison operation "<" valid in both directions.
-               IF ( inc * corr_index < inc * lc )  THEN 
-                  lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx &
-                       - coord_x(wall_index) ) / z0_l ) / logvelc1
-                  more = .TRUE.
-               ELSE
-                  lcr(alcorr) = 1.0_wp
-                  more = .FALSE.
-               ENDIF
-               lcorr  = lcorr + inc
-               alcorr = ABS( lcorr )
-            ENDDO
- 
-      END SELECT
-
-   END SUBROUTINE pmci_define_loglaw_correction_parameters
-
-
-
-   SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
-
-!
-!--   Finds the pivot node and te log-law factor for near-wall nodes for 
-!--   which the wall-parallel velocity components will be log-law corrected 
-!--   after interpolation. This subroutine is only for horizontal walls.
-!
-!--                               Antti Hellsten 30.12.2015
-      IMPLICIT NONE
-      REAL(wp),INTENT(OUT)  ::  logzc1     !:
-      REAL(wp), INTENT(IN)  ::  z0_l       !:
-      INTEGER(iwp), INTENT(IN)    ::  kb   !:
-      INTEGER(iwp), INTENT(OUT)   ::  lc   !:
+                IF ( corr_index < lc )  THEN
+                   lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
+                   more = .TRUE.
+                ELSE
+                   lcr(lcorr) = 1.0
+                   more = .FALSE.
+                ENDIF
+                lcorr = lcorr + 1
+             ENDDO
+
+          CASE (2)   !  j
+             more = .TRUE.
+             lcorr  = 0
+             alcorr = 0
+             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
+                corr_index = ij + lcorr   ! In this case (direction = 2) ij is j
+                IF ( lcorr == 0 )  THEN
+                   CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index,  &
+                                                z0_l, inc )
+                ENDIF
+
+!
+!--             The role of inc here is to make the comparison operation "<"
+!--             valid in both directions
+                IF ( inc * corr_index < inc * lc )  THEN
+                   lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy   &
+                                         - coord_y(wall_index) ) / z0_l )      &
+                                 / logvelc1
+                   more = .TRUE.
+                ELSE
+                   lcr(alcorr) = 1.0_wp
+                   more = .FALSE.
+                ENDIF
+                lcorr  = lcorr + inc
+                alcorr = ABS( lcorr )
+             ENDDO
+
+          CASE (3)   !  i
+             more = .TRUE.
+             lcorr  = 0
+             alcorr = 0
+             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
+                corr_index = ij + lcorr   ! In this case (direction = 3) ij is i
+                IF ( lcorr == 0 )  THEN
+                   CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index,  &
+                                                z0_l, inc )
+                ENDIF
+!
+!--             The role of inc here is to make the comparison operation "<"
+!--             valid in both directions
+                IF ( inc * corr_index < inc * lc )  THEN
+                   lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx   &
+                                         - coord_x(wall_index) ) / z0_l )      &
+                                 / logvelc1
+                   more = .TRUE.
+                ELSE
+                   lcr(alcorr) = 1.0_wp
+                   more = .FALSE.
+                ENDIF
+                lcorr  = lcorr + inc
+                alcorr = ABS( lcorr )
+             ENDDO
+
+       END SELECT
+
+    END SUBROUTINE pmci_define_loglaw_correction_parameters
+
+
+
+    SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
+!
+!--    Finds the pivot node and te log-law factor for near-wall nodes for
+!--    which the wall-parallel velocity components will be log-law corrected
+!--    after interpolation. This subroutine is only for horizontal walls.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN)  ::  kb   !:
+       INTEGER(iwp), INTENT(OUT) ::  lc   !:
+
+       INTEGER(iwp) ::  kbc    !:
+       INTEGER(iwp) ::  k1     !:
+
+       REAL(wp),INTENT(OUT) ::  logzc1     !:
+       REAL(wp), INTENT(IN) ::  z0_l       !:
+
+       REAL(wp) ::  zuc1   !:
+
+
+       kbc = nzb + 1
+!
+!--    kbc is the first coarse-grid point above the surface
+       DO  WHILE ( cg%zu(kbc) < zu(kb) )
+          kbc = kbc + 1
+       ENDDO
+       zuc1  = cg%zu(kbc)
+       k1    = kb + 1
+       DO  WHILE ( zu(k1) < zuc1 )  !  Important: must be <, not <=
+          k1 = k1 + 1
+       ENDDO
+       logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
+       lc = k1
+
+    END SUBROUTINE pmci_find_logc_pivot_k
+
+
+
+    SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
+!
+!--    Finds the pivot node and te log-law factor for near-wall nodes for
+!--    which the wall-parallel velocity components will be log-law corrected
+!--    after interpolation. This subroutine is only for vertical walls on
+!--    south/north sides of the node.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN)  ::  inc    !:  increment must be 1 or -1.
+       INTEGER(iwp), INTENT(IN)  ::  j      !:
+       INTEGER(iwp), INTENT(IN)  ::  jw     !:
+       INTEGER(iwp), INTENT(OUT) ::  lc     !:
+
+       INTEGER(iwp) ::  j1       !:
+
+       REAL(wp), INTENT(IN) ::  z0_l   !:
+
+       REAL(wp) ::  logyc1   !:
+       REAL(wp) ::  yc1      !:
+
+!
+!--    yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from
+!--    the wall
+       yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy
+!
+!--    j1 is the first fine-grid index further away from the wall than yc1
+       j1 = j
+!
+!--    Important: must be <, not <=
+       DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )
+          j1 = j1 + inc
+       ENDDO
+
+       logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
+       lc = j1
+
+    END SUBROUTINE pmci_find_logc_pivot_j
+
+
+
+    SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
+!
+!--    Finds the pivot node and the log-law factor for near-wall nodes for
+!--    which the wall-parallel velocity components will be log-law corrected
+!--    after interpolation. This subroutine is only for vertical walls on
+!--    south/north sides of the node.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp), INTENT(IN)  ::  i      !:
+       INTEGER(iwp), INTENT(IN)  ::  inc    !: increment must be 1 or -1.
+       INTEGER(iwp), INTENT(IN)  ::  iw     !:
+       INTEGER(iwp), INTENT(OUT) ::  lc     !:
+
+       INTEGER(iwp) ::  i1       !:
+
+       REAL(wp), INTENT(IN) ::  z0_l   !:
+
+       REAL(wp) ::  logxc1   !:
+       REAL(wp) ::  xc1      !:
+
+!
+!--    xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from
+!--    the wall
+       xc1  = coord_x(iw) + 0.5_wp * inc * cg%dx
+!
+!--    i1 is the first fine-grid index futher away from the wall than xc1.
+       i1 = i
+!
+!--    Important: must be <, not <=
+       DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc * xc1 )
+          i1 = i1 + inc
+       ENDDO
       
-      REAL(wp)     ::  zuc1   !:
-      INTEGER(iwp) ::  kbc    !:
-      INTEGER(iwp) ::  k1     !:
-
-
-      kbc = nzb + 1
-      DO  WHILE ( cg%zu(kbc) < zu(kb) )   !  kbc is the first coarse-grid point above the surface.
-         kbc = kbc + 1
-      ENDDO
-      zuc1  = cg%zu(kbc)
-      k1    = kb + 1  
-      DO  WHILE ( zu(k1) < zuc1 )         !  Important: must be <, not <=
-         k1 = k1 + 1
-      ENDDO
-      logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
-      lc = k1
-
-   END SUBROUTINE pmci_find_logc_pivot_k
-
-
-
-   SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
-
-!
-!--   Finds the pivot node and te log-law factor for near-wall nodes for 
-!--   which the wall-parallel velocity components will be log-law corrected 
-!--   after interpolation. This subroutine is only for vertical walls on 
-!--   south/north sides of the node.
-!
-!--                               Antti Hellsten 5.1.2016
-      IMPLICIT NONE
-      REAL(wp), INTENT(IN)      ::  z0_l   !:
-      INTEGER(iwp), INTENT(IN)  ::  inc    !:  increment must be 1 or -1.
-      INTEGER(iwp), INTENT(IN)  ::  j      !:
-      INTEGER(iwp), INTENT(IN)  ::  jw     !:
-      INTEGER(iwp), INTENT(OUT) ::  lc     !:
-      
-      REAL(wp)     ::  logyc1   !:
-      REAL(wp)     ::  yc1      !:
-      INTEGER(iwp) ::  j1       !:
-
-!
-!--   yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from the wall.
-      yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy  
-
-!
-!--   j1 is the first fine-grid index futher away from the wall than yc1.
-      j1 = j   
-      DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )   !  Important: must be <, not <=
-         j1 = j1 + inc
-      ENDDO
-
-      logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
-      lc = j1
-
-   END SUBROUTINE pmci_find_logc_pivot_j
-
-
-
-   SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
-
-!
-!--   Finds the pivot node and the log-law factor for near-wall nodes for 
-!--   which the wall-parallel velocity components will be log-law corrected 
-!--   after interpolation. This subroutine is only for vertical walls on 
-!--   south/north sides of the node.
-!
-!--                               Antti Hellsten 8.1.2016
-
-      IMPLICIT NONE
-      REAL(wp), INTENT(IN)      ::  z0_l   !:
-      INTEGER(iwp), INTENT(IN)  ::  i      !:
-      INTEGER(iwp), INTENT(IN)  ::  inc    !: increment must be 1 or -1.
-      INTEGER(iwp), INTENT(IN)  ::  iw     !:
-      INTEGER(iwp), INTENT(OUT) ::  lc     !:
-
-      REAL(wp)     ::  logxc1   !:
-      REAL(wp)     ::  xc1      !:
-      INTEGER(iwp) ::  i1       !:
-
-!
-!--   xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from the wall.
-      xc1  = coord_x(iw) + 0.5_wp *inc * cg%dx  
-
-!
-!--   i1 is the first fine-grid index futher away from the wall than xc1.
-      i1 = i   
-      DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc *xc1 )   !  Important: must be <, not <=
-         i1 = i1 + inc
-      ENDDO
-      
-      logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
-      lc = i1
-
-   END SUBROUTINE pmci_find_logc_pivot_i
-
-
-
+       logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
+       lc = i1
+
+    END SUBROUTINE pmci_find_logc_pivot_i
+
+
+!-- TO_DO:  indentation and wrap long lines from here on to the end of the file
    SUBROUTINE pmci_init_anterp_tophat
 !
@@ -1831,20 +1904,21 @@
 !--   corresponding coarse-grid array index and the 
 !--   Under-relaxation coefficients to be used by anterp_tophat.
-!
-!--                               Antti Hellsten 9.10.2015.
+
       IMPLICIT NONE
-      INTEGER(iwp) ::  i        !:
+
+      INTEGER(iwp) ::  i        !: Fine-grid index
+      INTEGER(iwp) ::  ii       !: Coarse-grid index
       INTEGER(iwp) ::  istart   !:
-      INTEGER(iwp) ::  j        !:
+      INTEGER(iwp) ::  j        !: Fine-grid index
+      INTEGER(iwp) ::  jj       !: Coarse-grid index
       INTEGER(iwp) ::  jstart   !:
-      INTEGER(iwp) ::  k        !:
+      INTEGER(iwp) ::  k        !: Fine-grid index
+      INTEGER(iwp) ::  kk       !: Coarse-grid index
       INTEGER(iwp) ::  kstart   !:
-      INTEGER(iwp) ::  l        !:
-      INTEGER(iwp) ::  m        !:
-      INTEGER(iwp) ::  n        !:
       REAL(wp)     ::  xi       !:
       REAL(wp)     ::  eta      !:
       REAL(wp)     ::  zeta     !:
      
+
 !
 !--   Default values:
@@ -1866,16 +1940,16 @@
 
 !
-!--   First determine kceu and kcew that are the coarse-grid upper bounds for index k.
-      n = 0
-      DO  WHILE ( cg%zu(n) < zu(nzt) )
-         n = n + 1
+!--   First determine kctu and kctw that are the coarse-grid upper bounds for index k.
+      kk = 0
+      DO  WHILE ( cg%zu(kk) < zu(nzt) )
+         kk = kk + 1
       ENDDO
-      kceu = n - 1 
-
-      n = 0
-      DO  WHILE ( cg%zw(n) < zw(nzt-1) )
-         n = n + 1
+      kctu = kk - 1 
+
+      kk = 0
+      DO  WHILE ( cg%zw(kk) < zw(nzt-1) )
+         kk = kk + 1
       ENDDO
-      kcew = n - 1 
+      kctw = kk - 1 
 
       ALLOCATE( iflu(icl:icr) )
@@ -1887,23 +1961,23 @@
       ALLOCATE( jfuv(jcs:jcn) )
       ALLOCATE( jfuo(jcs:jcn) )
-      ALLOCATE( kflw(0:kcew) )
-      ALLOCATE( kflo(0:kceu) )
-      ALLOCATE( kfuw(0:kcew) )
-      ALLOCATE( kfuo(0:kceu) )
+      ALLOCATE( kflw(0:kctw) )
+      ALLOCATE( kflo(0:kctu) )
+      ALLOCATE( kfuw(0:kctw) )
+      ALLOCATE( kfuo(0:kctu) )
 
 !
 !--   i-indices of u for each l-index value.    
       istart = nxlg
-      DO  l = icl, icr  
+      DO  ii = icl, icr  
          i = istart
-         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
+         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
             i = i + 1
          ENDDO
-         iflu(l) = MIN( MAX( i, nxlg ), nxrg )
-         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
+         iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
+         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
             i = i + 1
          ENDDO
-         ifuu(l) = MIN( MAX( i, nxlg ), nxrg )
-         istart = iflu(l)
+         ifuu(ii) = MIN( MAX( i, nxlg ), nxrg )
+         istart = iflu(ii)
       ENDDO
 
@@ -1911,15 +1985,15 @@
 !--   i-indices of others for each l-index value.
       istart = nxlg
-      DO  l = icl, icr   
+      DO  ii = icl, icr   
          i = istart
-         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) )  .AND.  ( i < nxrg ) ) 
+         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(ii) )  .AND.  ( i < nxrg ) ) 
             i = i + 1
          ENDDO
-         iflo(l) = MIN( MAX( i, nxlg ), nxrg )
-         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) + cg%dx )  .AND.  ( i < nxrg ) ) 
+         iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
+         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(ii) + cg%dx )  .AND.  ( i < nxrg ) ) 
             i = i + 1
          ENDDO
-         ifuo(l) = MIN(MAX(i,nxlg),nxrg)
-         istart = iflo(l)
+         ifuo(ii) = MIN(MAX(i,nxlg),nxrg)
+         istart = iflo(ii)
       ENDDO
 
@@ -1927,15 +2001,15 @@
 !--   j-indices of v for each m-index value.
       jstart = nysg
-      DO  m = jcs, jcn
+      DO  jj = jcs, jcn
          j = jstart
-         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
+         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
             j = j + 1
          ENDDO
-         jflv(m) = MIN( MAX( j, nysg ), nyng ) 
-         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
+         jflv(jj) = MIN( MAX( j, nysg ), nyng ) 
+         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
             j = j + 1
          ENDDO
-         jfuv(m) = MIN( MAX( j, nysg ), nyng )
-         jstart = jflv(m)
+         jfuv(jj) = MIN( MAX( j, nysg ), nyng )
+         jstart = jflv(jj)
       ENDDO
 
@@ -1943,15 +2017,15 @@
 !--   j-indices of others for each m-index value.
       jstart = nysg
-      DO  m = jcs, jcn
+      DO  jj = jcs, jcn
          j = jstart
-         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) )  .AND.  ( j < nyng ) )  
+         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(jj) )  .AND.  ( j < nyng ) )  
             j = j + 1
          ENDDO
-         jflo(m) = MIN( MAX( j, nysg ), nyng )
-         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) + cg%dy )  .AND.  ( j < nyng ) ) 
+         jflo(jj) = MIN( MAX( j, nysg ), nyng )
+         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(jj) + cg%dy )  .AND.  ( j < nyng ) ) 
             j = j + 1
          ENDDO
-         jfuo(m) = MIN( MAX( j, nysg ), nyng )
-         jstart = jflv(m)
+         jfuo(jj) = MIN( MAX( j, nysg ), nyng )
+         jstart = jflv(jj)
       ENDDO
 
@@ -1961,15 +2035,15 @@
       kflw(0) = 0
       kfuw(0) = 0
-      DO  n = 1, kcew
+      DO  kk = 1, kctw
          k = kstart
-         DO  WHILE ( ( zw(k) < cg%zw(n) - 0.5_wp * cg%dzw(n) )  .AND.  ( k < nzt ) ) 
+         DO  WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) )  .AND.  ( k < nzt ) ) 
             k = k + 1
          ENDDO
-         kflw(n) = MIN( MAX( k, 1 ), nzt + 1 )
-         DO  WHILE ( ( zw(k) < cg%zw(n) + 0.5_wp * cg%dzw(n+1) )  .AND.  ( k < nzt ) ) 
+         kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+         DO  WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) )  .AND.  ( k < nzt ) ) 
             k = k + 1
          ENDDO
-         kfuw(n) = MIN( MAX( k, 1 ), nzt + 1 )
-         kstart = kflw(n)
+         kfuw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+         kstart = kflw(kk)
       ENDDO
 
@@ -1979,15 +2053,15 @@
       kflo(0) = 0
       kfuo(0) = 0
-      DO  n = 1, kceu
+      DO  kk = 1, kctu
          k = kstart
-         DO  WHILE ( ( zu(k) < cg%zu(n) - 0.5_wp * cg%dzu(n) )  .AND.  ( k < nzt ) ) 
+         DO  WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) )  .AND.  ( k < nzt ) ) 
             k = k + 1
          ENDDO
-         kflo(n) = MIN( MAX( k, 1 ), nzt + 1 )
-         DO  WHILE ( ( zu(k) < cg%zu(n) + 0.5_wp * cg%dzu(n+1) )  .AND.  ( k < nzt ) )  
+         kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+         DO  WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) )  .AND.  ( k < nzt ) )  
             k = k + 1
          ENDDO
-         kfuo(n) = MIN( MAX( k-1, 1 ), nzt + 1 )
-         kstart = kflo(n)
+         kfuo(kk) = MIN( MAX( k-1, 1 ), nzt + 1 )
+         kstart = kflo(kk)
       ENDDO
      
@@ -1995,31 +2069,33 @@
 !--   Spatial under-relaxation coefficients
       ALLOCATE( frax(icl:icr) )
-      DO  l = icl, icr
+
+      DO  ii = icl, icr
          IF ( nest_bound_l ) THEN
-            xi   = ( ( cg%coord_x(l) - lower_left_coord_x ) / anterp_relax_length_l )**4
+            xi   = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) / anterp_relax_length_l )**4
          ELSEIF ( nest_bound_r ) THEN
-            xi   = ( ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(l) ) / anterp_relax_length_r )**4
+            xi   = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(ii) ) ) / anterp_relax_length_r )**4
          ELSE
             xi   = 999999.9_wp
          ENDIF
-         frax(l)  = xi / ( 1.0_wp + xi )
+         frax(ii)  = xi / ( 1.0_wp + xi )
       ENDDO
 
       ALLOCATE( fray(jcs:jcn) )
-      DO  m = jcs, jcn
+
+      DO  jj = jcs, jcn
          IF ( nest_bound_s ) THEN           
-            eta  = ( ( cg%coord_y(m) - lower_left_coord_y ) / anterp_relax_length_s )**4
+            eta  = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) / anterp_relax_length_s )**4
          ELSEIF ( nest_bound_n ) THEN
-            eta  = ( (lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(m)) / anterp_relax_length_n )**4
+            eta  = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(jj)) ) / anterp_relax_length_n )**4
          ELSE
             eta  = 999999.9_wp
          ENDIF
-         fray(m)  = eta / ( 1.0_wp + eta )
+         fray(jj)  = eta / ( 1.0_wp + eta )
       ENDDO
       
-      ALLOCATE( fraz(0:kceu) )
-      DO  n = 0, kceu        
-         zeta = ( ( zu(nzt) - cg%zu(n) ) / anterp_relax_length_t )**4       
-         fraz(n) = zeta / ( 1.0_wp + zeta )
+      ALLOCATE( fraz(0:kctu) )
+      DO  kk = 0, kctu        
+         zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4       
+         fraz(kk) = zeta / ( 1.0_wp + zeta )
       ENDDO
 
@@ -2173,23 +2249,544 @@
 
 
+
+ SUBROUTINE pmci_set_array_pointer( name, client_id, nz_cl )
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)          ::  client_id   !:
+    INTEGER, INTENT(IN)          ::  nz_cl       !:
+    CHARACTER(LEN=*), INTENT(IN) ::  name        !:
+
+#if defined( __parallel )
+    INTEGER(iwp) ::  ierr        !:
+    INTEGER(iwp) ::  istat       !:
+
+    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d        !:
+    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d_sec    !:
+    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d        !:
+    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d_sec    !:
+
+
+    NULLIFY( p_3d )
+    NULLIFY( p_2d )
+
+!
+!-- List of array names, which can be coupled.
+!-- In case of 3D please change also the second array for the pointer version
+    IF ( TRIM(name) == "u"  )  p_3d => u
+    IF ( TRIM(name) == "v"  )  p_3d => v
+    IF ( TRIM(name) == "w"  )  p_3d => w
+    IF ( TRIM(name) == "e"  )  p_3d => e
+    IF ( TRIM(name) == "pt" )  p_3d => pt
+    IF ( TRIM(name) == "q"  )  p_3d => q
+!
+!-- 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 defined( __nopointer )
+    IF ( ASSOCIATED( p_3d ) )  THEN
+       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz )
+    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
+       CALL pmc_s_set_dataarray( client_id, p_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'
+          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
+       ELSE
+!
+!--       Avoid others to continue
+          CALL MPI_BARRIER( comm2d, ierr )
+       ENDIF
+    ENDIF
+#else
+    IF ( TRIM(name) == "u"  )  p_3d_sec => u_2
+    IF ( TRIM(name) == "v"  )  p_3d_sec => v_2
+    IF ( TRIM(name) == "w"  )  p_3d_sec => w_2
+    IF ( TRIM(name) == "e"  )  p_3d_sec => e_2
+    IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
+    IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
+
+    IF ( ASSOCIATED( p_3d ) )  THEN
+       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz, &
+                                 array_2 = p_3d_sec )
+    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
+       CALL pmc_s_set_dataarray( client_id, p_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'
+          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
+       ELSE
+!
+!--       Avoid others to continue
+          CALL MPI_BARRIER( comm2d, ierr )
+       ENDIF
+
+   ENDIF
+#endif
+
+#endif
+ END SUBROUTINE pmci_set_array_pointer
+
+
+
+ SUBROUTINE pmci_create_client_arrays( name, is, ie, js, je, nzc  )
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN) ::  name    !:
+
+    INTEGER(iwp), INTENT(IN) ::  ie      !:
+    INTEGER(iwp), INTENT(IN) ::  is      !:
+    INTEGER(iwp), INTENT(IN) ::  je      !:
+    INTEGER(iwp), INTENT(IN) ::  js      !:
+    INTEGER(iwp), INTENT(IN) ::  nzc     !:  Note that nzc is cg%nz
+
+#if defined( __parallel )
+    INTEGER(iwp) ::  ierr    !:
+    INTEGER(iwp) ::  istat   !:
+
+    REAL(wp), POINTER,DIMENSION(:,:)   ::  p_2d    !:
+    REAL(wp), POINTER,DIMENSION(:,:,:) ::  p_3d    !:
+
+
+    NULLIFY( p_3d )
+    NULLIFY( p_2d )
+
+!
+!-- List of array names, which can be coupled.
+!-- TO_DO: Antti: what is the meaning of the next line?
+!-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1.
+    IF ( TRIM( name ) == "u" )  THEN
+       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
+       p_3d => uc
+    ELSEIF ( TRIM( name ) == "v" )  THEN
+       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1, js:je, is:ie) )
+       p_3d => vc
+    ELSEIF ( TRIM( name ) == "w" )  THEN
+       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1, js:je, is:ie) )
+       p_3d => wc
+    ELSEIF ( TRIM( name ) == "e" )  THEN
+       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1, js:je, is:ie) )
+       p_3d => ec
+    ELSEIF ( TRIM( name ) == "pt")  THEN
+       IF ( .NOT. ALLOCATED( ptc ) ) ALLOCATE( ptc(0:nzc+1, js:je, is:ie) )
+       p_3d => ptc
+    ELSEIF ( TRIM( name ) == "q")  THEN
+       IF ( .NOT. ALLOCATED( qc ) ) ALLOCATE( qc(0:nzc+1, js:je, is:ie) )
+       p_3d => qc
+    !ELSEIF (trim(name) == "z0") then
+       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
+       !p_2d => z0c
+    ENDIF
+
+    IF ( ASSOCIATED( p_3d ) )  THEN
+       CALL pmc_c_set_dataarray( p_3d )
+    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
+       CALL pmc_c_set_dataarray( p_2d )
+    ELSE
+!
+!--    Give only one message for the first client domain
+       IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
+
+          message_string = 'pointer for array "' // TRIM( name ) //            &
+                           '" can''t be associated'
+          CALL message( 'pmci_create_client_arrays', 'PA0170', 3, 2, 0, 6, 0 )
+       ELSE
+!
+!--       Avoid others to continue
+          CALL MPI_BARRIER( comm2d, ierr )
+       ENDIF
+    ENDIF
+
+#endif
+ END SUBROUTINE pmci_create_client_arrays
+
+
+
+ SUBROUTINE pmci_server_initialize
+!-- TO_DO: add general explanations about what this subroutine does
+#if defined( __parallel )
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  client_id   !:
+    INTEGER(iwp) ::  m           !:
+
+    REAL(wp) ::  waittime    !:
+
+
+    DO  m = 1, SIZE( pmc_server_for_client ) - 1
+       client_id = pmc_server_for_client(m)
+       CALL pmc_s_fillbuffer( client_id, waittime=waittime )
+    ENDDO
+
+#endif
+ END SUBROUTINE pmci_server_initialize
+
+
+
+ SUBROUTINE pmci_client_initialize
+!-- TO_DO: add general explanations about what this subroutine does
+#if defined( __parallel )
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i          !:
+    INTEGER(iwp) ::  icl        !:
+    INTEGER(iwp) ::  icr        !:
+    INTEGER(iwp) ::  j          !:
+    INTEGER(iwp) ::  jcn        !:
+    INTEGER(iwp) ::  jcs        !:
+
+    REAL(wp) ::  waittime   !:
+
+!
+!-- Root id is never a client
+    IF ( cpl_id > 1 )  THEN
+
+!
+!--    Client domain boundaries in the server index space
+       icl = coarse_bound(1)
+       icr = coarse_bound(2)
+       jcs = coarse_bound(3)
+       jcn = coarse_bound(4)
+
+!
+!--    Get data from server
+       CALL pmc_c_getbuffer( waittime = waittime )
+
+!
+!--    The interpolation.
+       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
+                                   r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
+       CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,   &
+                                   r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
+       CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,   &
+                                   r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
+       CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
+                                   r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
+       CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
+                                   r2yo, r1zo, r2zo, nzb_s_inner, 's' )
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          CALL pmci_interp_tril_all ( q, qc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                      r2yo, r1zo, r2zo, nzb_s_inner, 's' )
+       ENDIF
+
+       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
+                u(nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
+                v(nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
+                w(nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
+                e(nzb:nzb_s_inner(j,i),j,i)   = 0.0_wp
+                u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
+                v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
+                w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
+                e_p(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
+
+
+ CONTAINS
+
+
+    SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
+                                     r1z, r2z, kb, var )
+!
+!--    Interpolation of the internal values for the client-domain initialization
+!--    This subroutine is based on trilinear interpolation.
+!--    Coding based on interp_tril_lr/sn/t
+       IMPLICIT NONE
+
+       CHARACTER(LEN=1), INTENT(IN) :: var  !:
+
+       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !:
+       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !:
+       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
+       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !:
+
+       INTEGER(iwp) ::  i      !:
+       INTEGER(iwp) ::  ib     !:
+       INTEGER(iwp) ::  ie     !:
+       INTEGER(iwp) ::  j      !:
+       INTEGER(iwp) ::  jb     !:
+       INTEGER(iwp) ::  je     !:
+       INTEGER(iwp) ::  k      !:
+       INTEGER(iwp) ::  k1     !:
+       INTEGER(iwp) ::  kbc    !:
+       INTEGER(iwp) ::  l      !:
+       INTEGER(iwp) ::  m      !:
+       INTEGER(iwp) ::  n      !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
+       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc    !:
+       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x   !:
+       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x   !:
+       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y   !:
+       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y   !:
+       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z   !:
+       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z   !:
+
+       REAL(wp) ::  fk         !:
+       REAL(wp) ::  fkj        !:
+       REAL(wp) ::  fkjp       !:
+       REAL(wp) ::  fkp        !:
+       REAL(wp) ::  fkpj       !:
+       REAL(wp) ::  fkpjp      !:
+       REAL(wp) ::  logratio   !:
+       REAL(wp) ::  logzuc1    !:
+       REAL(wp) ::  zuc1       !:
+
+
+       ib = nxl
+       ie = nxr
+       jb = nys
+       je = nyn
+       IF ( nest_bound_l )  THEN
+          ib = nxl - 1
+!
+!--       For u, nxl is a ghost node, but not for the other variables
+          IF ( var == 'u' )  THEN
+             ib = nxl
+          ENDIF
+       ENDIF
+       IF ( nest_bound_s )  THEN
+          jb = nys - 1
+!
+!--       For v, nys is a ghost node, but not for the other variables
+          IF ( var == 'v' )  THEN
+             jb = nys
+          ENDIF
+       ENDIF
+       IF ( nest_bound_r )  THEN
+          ie = nxr + 1
+       ENDIF
+       IF ( nest_bound_n )  THEN
+          je = nyn + 1
+       ENDIF
+
+!
+!--    Trilinear interpolation.
+       DO  i = ib, ie
+          DO  j = jb, je
+             DO  k = kb(j,i), nzt + 1
+                l = ic(i)
+                m = jc(j)
+                n = kc(k)
+                fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
+                fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
+                fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
+                fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
+                fk       = r1y(j) * fkj  + r2y(j) * fkjp
+                fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
+                f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Correct the interpolated values of u and v in near-wall nodes, i.e. in
+!--    the nodes below the coarse-grid nodes with k=1. The corrction is only
+!--    made over horizontal wall surfaces in this phase. For the nest boundary
+!--    conditions, a corresponding correction is made for all vertical walls,
+!--    too.
+       IF ( var == 'u' .OR. var == 'v' )  THEN
+          DO  i = ib, nxr
+             DO  j = jb, nyn
+                kbc = 1
+!
+!--             kbc is the first coarse-grid point above the surface
+                DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )
+                   kbc = kbc + 1
+                ENDDO
+                zuc1 = cg%zu(kbc)
+                k1   = kb(j,i) + 1
+                DO  WHILE ( zu(k1) < zuc1 )
+                   k1 = k1 + 1
+                ENDDO
+                logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
+
+                k = kb(j,i) + 1
+                DO  WHILE ( zu(k) < zuc1 )
+                   logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) / logzuc1
+                   f(k,j,i) = logratio * f(k1,j,i)
+                   k  = k + 1
+                ENDDO
+                f(kb(j,i),j,i) = 0.0_wp
+             ENDDO
+          ENDDO
+
+       ELSEIF ( var == 'w' )  THEN
+
+          DO  i = ib, nxr
+              DO  j = jb, nyn
+                f(kb(j,i),j,i) = 0.0_wp
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+    END SUBROUTINE pmci_interp_tril_all
+
+#endif
+ END SUBROUTINE pmci_client_initialize
+
+
+
+ SUBROUTINE pmci_ensure_nest_mass_conservation
+
+#if defined( __parallel )
+!
+!-- Adjust the volume-flow rate through the top boundary so that the net volume
+!-- flow through all boundaries of the current nest domain becomes zero.
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  i                          !:
+    INTEGER(iwp) ::  ierr                       !:
+    INTEGER(iwp) ::  j                          !:
+    INTEGER(iwp) ::  k                          !:
+
+    REAL(wp) ::  dxdy                            !:
+    REAL(wp) ::  innor                           !:
+    REAL(wp) ::  w_lt                            !:
+    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !:
+
+!
+!-- Sum up the volume flow through the left/right boundaries
+    volume_flow(1)   = 0.0_wp
+    volume_flow_l(1) = 0.0_wp
+
+    IF ( nest_bound_l )  THEN
+       i = 0
+       innor = dy
+       DO   j = nys, nyn
+          DO   k = nzb_u_inner(j,i)+1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
+          ENDDO
+       ENDDO
+    ENDIF
+
+    IF ( nest_bound_r )  THEN
+       i = nx + 1
+       innor = -dy
+       DO   j = nys, nyn
+          DO   k = nzb_u_inner(j,i)+1, nzt
+             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
+          ENDDO
+       ENDDO
+    ENDIF
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
+                        MPI_SUM, comm2d, ierr )
+#else
+    volume_flow(1) = volume_flow_l(1)
+#endif
+
+!
+!-- Sum up the volume flow through the south/north boundaries
+    volume_flow(2)   = 0.0_wp
+    volume_flow_l(2) = 0.0_wp
+
+    IF ( nest_bound_s )  THEN
+       j = 0
+       innor = dx
+       DO   i = nxl, nxr
+          DO   k = nzb_v_inner(j,i)+1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
+          ENDDO
+       ENDDO
+    ENDIF
+
+    IF ( nest_bound_n )  THEN
+       j = ny + 1
+       innor = -dx
+       DO   i = nxl, nxr
+          DO   k = nzb_v_inner(j,i)+1, nzt
+             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
+          ENDDO
+       ENDDO
+    ENDIF
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL,         &
+                        MPI_SUM, comm2d, ierr )
+#else
+    volume_flow(2) = volume_flow_l(2)
+#endif
+
+!
+!-- Sum up the volume flow through the top boundary
+    volume_flow(3)   = 0.0_wp
+    volume_flow_l(3) = 0.0_wp
+    dxdy = dx * dy
+    k = nzt
+    DO   i = nxl, nxr
+       DO   j = nys, nyn
+          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
+       ENDDO
+    ENDDO
+
+#if defined( __parallel )
+    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL,         &
+                        MPI_SUM, comm2d, ierr )
+#else
+    volume_flow(3) = volume_flow_l(3)
+#endif
+
+!
+!-- Correct the top-boundary value of w
+    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
+    DO   i = nxl, nxr
+       DO   j = nys, nyn
+          DO  k = nzt, nzt + 1
+             w(k,j,i) = w(k,j,i) + w_lt
+          ENDDO
+       ENDDO
+    ENDDO
+
+#endif
+ END SUBROUTINE pmci_ensure_nest_mass_conservation
+
+
+!-- TO_DO: the timestep sycnchronization could be done easier using
+!--        an MPI_ALLREDUCE with MIN over MPI_COMM_WORLD
  SUBROUTINE pmci_server_synchronize
 
 #if defined( __parallel )
 !
-!-- Unify the time steps for each model and synchronize.
-!-- This is based on the assumption that the native time step
-!-- (original dt_3d) of any server is always larger than the smallest 
-!-- native time step of it s clients. 
+!-- Unify the time steps for each model and synchronize. This is based on the
+!-- assumption that the native time step (original dt_3d) of any server is
+!-- always larger than the smallest native time step of it s clients.
     IMPLICIT NONE
-    INTEGER(iwp)           ::  client_id   !:
+
+    INTEGER(iwp) ::  client_id   !:
+    INTEGER(iwp) ::  ierr        !:
+    INTEGER(iwp) ::  m           !:
+
     REAL(wp), DIMENSION(1) ::  dtc         !:
     REAL(wp), DIMENSION(1) ::  dtl         !:
-    INTEGER(iwp)           ::  ierr        !:
-    INTEGER(iwp)           ::  m           !:    
-
-!
-!-- First find the smallest native time step of all the clients of the current server.
+
+!
+!-- First find the smallest native time step of all the clients of the current
+!-- server.
     dtl(1) = 999999.9_wp
-    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
+    DO  m = 1, SIZE( PMC_Server_for_Client )-1
        client_id = PMC_Server_for_Client(m)
        IF ( myid == 0 )  THEN
@@ -2201,9 +2798,9 @@
 
 !
-!-- Broadcast the unified time step to all server processes.
+!-- Broadcast the unified time step to all server processes
     CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
 
 !
-!-- Send the new time step to all the clients of the current server.
+!-- Send the new time step to all the clients of the current server
     DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
        client_id = PMC_Server_for_Client(m)
@@ -2222,17 +2819,18 @@
 #if defined( __parallel )
 !
-!-- Unify the time steps for each model and synchronize.
-!-- This is based on the assumption that the native time step
-!-- (original dt_3d) of any server is always larger than the smallest 
-!-- native time step of it s clients. 
+!-- Unify the time steps for each model and synchronize. This is based on the
+!-- assumption that the native time step (original dt_3d) of any server is
+!-- always larger than the smallest native time step of it s clients.
 
     IMPLICIT NONE
+
+    INTEGER(iwp) ::  ierr   !:
+
     REAL(wp), DIMENSION(1) ::  dtl    !:
     REAL(wp), DIMENSION(1) ::  dts    !:
-    INTEGER(iwp)           ::  ierr   !:
     
 
     dtl(1) = dt_3d
-    IF ( cpl_id > 1 )  THEN   !  Root id is never a client
+    IF ( cpl_id > 1 )  THEN
        IF ( myid==0 )  THEN
           CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
@@ -2242,5 +2840,5 @@
 
 !
-!--    Broadcast the unified time step to all server processes.
+!--    Broadcast the unified time step to all server processes
        CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
     ENDIF
@@ -2255,5 +2853,6 @@
     IMPLICIT NONE
 
-    INTEGER(iwp),INTENT(IN) ::  swaplevel  !: swaplevel (1 or 2) of PALM's timestep
+    INTEGER(iwp),INTENT(IN) ::  swaplevel  !: swaplevel (1 or 2) of PALM's
+                                           !: timestep
 
     INTEGER(iwp)            ::  client_id  !:
@@ -2278,18 +2877,20 @@
 
 #if defined( __parallel )
-    INTEGER(iwp)            ::  client_id   !:
-    INTEGER(iwp)            ::  i           !:
-    INTEGER(iwp)            ::  j           !:
-    INTEGER(iwp)            ::  ierr        !:
-    INTEGER(iwp)            ::  m           !:
-    REAL(wp)                ::  waittime    !:
-    REAL(wp), DIMENSION(1)  ::  dtc         !:
-    REAL(wp), DIMENSION(1)  ::  dtl         !:
-
-!
-!-- First find the smallest native time step of all the clients of the current server.
+    INTEGER(iwp) ::  client_id   !:
+    INTEGER(iwp) ::  i           !:
+    INTEGER(iwp) ::  j           !:
+    INTEGER(iwp) ::  ierr        !:
+    INTEGER(iwp) ::  m           !:
+
+    REAL(wp)               ::  waittime    !:
+    REAL(wp), DIMENSION(1) ::  dtc         !:
+    REAL(wp), DIMENSION(1) ::  dtl         !:
+
+!
+!-- First find the smallest native time step of all the clients of the current
+!-- server.
     dtl(1) = 999999.9_wp
-    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
-       client_id = PMC_Server_for_Client(m)
+    DO  m = 1, SIZE( pmc_server_for_client )-1
+       client_id = pmc_server_for_client(m)
        IF ( myid==0 )  THEN
           CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )
@@ -2300,36 +2901,39 @@
 
 !
-!-- Broadcast the unified time step to all server processes.
+!-- Broadcast the unified time step to all server processes
     CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
 
-    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
+    DO  m = 1, SIZE( PMC_Server_for_Client )-1
        client_id = PMC_Server_for_Client(m)
-       CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
-
-!
-!--    Send the new time step to all the clients of the current server.
-       IF ( myid == 0 ) THEN
+       CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )
+
+!
+!--    Send the new time step to all the clients of the current server
+       IF ( myid == 0 )  THEN
           CALL pmc_send_to_client( client_id, dtl, SIZE( dtl ), 0, 102, ierr )
        ENDIF
-       CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
+       CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )
        
        IF ( direction == server_to_client )  THEN
-          CALL cpu_log( log_point_s(71), 'PMC Server Send', 'start' )
+          CALL cpu_log( log_point_s(71), 'pmc server send', 'start' )
           CALL pmc_s_fillbuffer( client_id, waittime=waittime )
-          CALL cpu_log( log_point_s(71), 'PMC Server Send', 'stop' )
-       ELSE   !  Communication from client to server.
-          CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'start' )
+          CALL cpu_log( log_point_s(71), 'pmc server send', 'stop' )
+       ELSE
+!
+!--       Communication from client to server
+          CALL cpu_log( log_point_s(72), 'pmc server recv', 'start' )
           client_id = pmc_server_for_client(m)
           CALL pmc_s_getdata_from_buffer( client_id )
-          CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'stop' )        
-
-!
-!--       The anterpolated data is now available in u etc.
+          CALL cpu_log( log_point_s(72), 'pmc server recv', 'stop' )
+
+!
+!--       The anterpolated data is now available in u etc
           IF ( topography /= 'flat' )  THEN
 
 !
 !--          Inside buildings/topography reset velocities and TKE back to zero.
-!--          Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
-!--          maybe revise later.
+!--          TO_DO: at least temperature should be included here immediately
+!--          Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at
+!--          present, maybe revise later.
              DO   i = nxlg, nxrg
                 DO   j = nysg, nyng
@@ -2337,5 +2941,5 @@
                    v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
                    w(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
-                   e(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
+                   e(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
                 ENDDO
              ENDDO
@@ -2356,18 +2960,19 @@
 
 #if defined( __parallel )
-    INTEGER(iwp)             ::  ierr        !:
-    INTEGER(iwp)             ::  icl         !:
-    INTEGER(iwp)             ::  icr         !:
-    INTEGER(iwp)             ::  jcs         !:
-    INTEGER(iwp)             ::  jcn         !:
+    INTEGER(iwp) ::  ierr        !:
+    INTEGER(iwp) ::  icl         !:
+    INTEGER(iwp) ::  icr         !:
+    INTEGER(iwp) ::  jcs         !:
+    INTEGER(iwp) ::  jcn         !:
     
-    REAL(wp), DIMENSION(1)   ::  dtl         !:
-    REAL(wp), DIMENSION(1)   ::  dts         !:
-    REAL(wp)                 ::  waittime    !:
+    REAL(wp) ::  waittime    !:
+
+    REAL(wp), DIMENSION(1) ::  dtl         !:
+    REAL(wp), DIMENSION(1) ::  dts         !:
 
 
     dtl = dt_3d
-    IF ( cpl_id > 1 )  THEN   !  Root id is never a client      
-       CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
+    IF ( cpl_id > 1 )  THEN
+       CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )
        IF ( myid==0 )  THEN
           CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
@@ -2379,5 +2984,5 @@
 !--    Broadcast the unified time step to all server processes.
        CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
-       CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
+       CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )
 
 !
@@ -2389,17 +2994,20 @@
 
        IF ( direction == server_to_client )  THEN
-          CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'start' )
+
+          CALL cpu_log( log_point_s(73), 'pmc client recv', 'start' )
           CALL pmc_c_getbuffer( WaitTime = WaitTime )
-          CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'stop' )
+          CALL cpu_log( log_point_s(73), 'pmc client recv', 'stop' )
 
           CALL pmci_interpolation
 
-       ELSE    !  IF ( direction == server_to_client )
-
+       ELSE
+!
+!--       direction == server_to_client
           CALL pmci_anterpolation
 
-          CALL cpu_log( log_point_s(74), 'PMC Client Send', 'start' )
+          CALL cpu_log( log_point_s(74), 'pmc client send', 'start' )
           CALL pmc_c_putbuffer( WaitTime = WaitTime )
-          CALL cpu_log( log_point_s(74), 'PMC Client Send', 'stop' )
+          CALL cpu_log( log_point_s(74), 'pmc client send', 'stop' )
+
        ENDIF
     ENDIF
@@ -2407,56 +3015,85 @@
  CONTAINS
 
-
-   SUBROUTINE pmci_interpolation
-
-!
-!--   A wrapper routine for all interpolation and extrapolation actions.
-      IMPLICIT NONE
-
-!
-!--   Add IF-condition here: IF not vertical nesting
-      IF ( nest_bound_l )  THEN   !  Left border pe
-         CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_u_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'u' )
-         CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
-                                   nzb_v_inner, logc_v_l, logc_ratio_v_l, nzt_topo_nestbc_l, 'l', 'v' )
-         CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
-                                   nzb_w_inner, logc_w_l, logc_ratio_w_l, nzt_topo_nestbc_l, 'l', 'w' )
-         CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'e' )
-         CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
-         IF ( humidity  .OR.  passive_scalar )  THEN
-            CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
-                                      nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
-         ENDIF
-         IF ( nesting_mode == 'one-way' )  THEN
-            CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
-            CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
-            CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
-            CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
-            CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
-            IF ( humidity  .OR.  passive_scalar )  THEN
-               CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
-            ENDIF
-         ENDIF
-      ENDIF
-      IF ( nest_bound_r )  THEN   !  Right border pe
-         CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_u_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'u' )
-         CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
-                                   nzb_v_inner, logc_v_r, logc_ratio_v_r, nzt_topo_nestbc_r, 'r', 'v' )
-         CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
-                                   nzb_w_inner, logc_w_r, logc_ratio_w_r, nzt_topo_nestbc_r, 'r', 'w' )
-         CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'e' )
-         CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                   nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
-         IF ( humidity  .OR.  passive_scalar )  THEN
-            CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
-                                      nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
-         ENDIF
-         IF ( nesting_mode == 'one-way' )  THEN
-            CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
+    SUBROUTINE pmci_interpolation
+
+!
+!--    A wrapper routine for all interpolation and extrapolation actions
+       IMPLICIT NONE
+
+!
+!--    Add IF-condition here: IF not vertical nesting
+!--    Left border pe:
+       IF ( nest_bound_l )  THEN
+          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_l,   &
+                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
+                                    'u' )
+          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
+                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_l,   &
+                                    logc_ratio_v_l, nzt_topo_nestbc_l, 'l',    &
+                                    'v' )
+          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_l,   &
+                                    logc_ratio_w_l, nzt_topo_nestbc_l, 'l',    &
+                                    'w' )
+          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_l,   &
+                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
+                                    'e' )
+          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_l,   &
+                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
+                                    's' )
+          IF ( humidity  .OR.  passive_scalar )  THEN
+             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
+                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
+                                       logc_u_l, logc_ratio_u_l,               &
+                                       nzt_topo_nestbc_l, 'l', 's' )
+          ENDIF
+
+          IF ( nesting_mode == 'one-way' )  THEN
+             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
+             CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
+             CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
+             CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
+             CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
+             IF ( humidity  .OR.  passive_scalar )  THEN
+                CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
+             ENDIF
+          ENDIF
+
+       ENDIF
+!
+!--    Right border pe
+       IF ( nest_bound_r )  THEN
+          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_r,   &
+                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
+                                    'u' )
+          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
+                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_r,   &
+                                    logc_ratio_v_r, nzt_topo_nestbc_r, 'r',    &
+                                    'v' )
+          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_r,   &
+                                    logc_ratio_w_r, nzt_topo_nestbc_r, 'r',    &
+                                    'w' )
+          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_r,   &
+                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
+                                    'e' )
+          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_r,   &
+                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
+                                    's' )
+          IF ( humidity  .OR.  passive_scalar )  THEN
+             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
+                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
+                                       logc_u_r, logc_ratio_u_r,               &
+                                       nzt_topo_nestbc_r, 'r', 's' )
+          ENDIF
+
+          IF ( nesting_mode == 'one-way' )  THEN
+             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
              CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'r', 'v' )
              CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'r', 'w' )
@@ -2467,20 +3104,36 @@
              ENDIF
           ENDIF
-       ENDIF
-       IF ( nest_bound_s )  THEN   !  South border pe
-          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
-                                    nzb_u_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'u' )
-          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
-                                    nzb_v_inner, logc_v_s, logc_ratio_v_s, nzt_topo_nestbc_s, 's', 'v' )
-          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
-                                    nzb_w_inner, logc_w_s, logc_ratio_w_s, nzt_topo_nestbc_s, 's', 'w' )
-          CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                    nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'e' )
-          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                    nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
+
+       ENDIF
+!
+!--    South border pe
+       IF ( nest_bound_s )  THEN
+          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_s,   &
+                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
+                                    'u' )
+          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
+                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_s,   &
+                                    logc_ratio_v_s, nzt_topo_nestbc_s, 's',    &
+                                    'v' )
+          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_s,   &
+                                    logc_ratio_w_s, nzt_topo_nestbc_s, 's',    &
+                                    'w' )
+          CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_s,   &
+                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
+                                    'e' )
+          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_s,   &
+                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
+                                    's' )
           IF ( humidity  .OR.  passive_scalar )  THEN
-             CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
-                                       nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
+             CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
+                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
+                                       logc_u_s, logc_ratio_u_s,               &
+                                       nzt_topo_nestbc_s, 's', 's' )
           ENDIF
+
           IF ( nesting_mode == 'one-way' )  THEN
              CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 's', 'u' )
@@ -2493,20 +3146,36 @@
              ENDIF
           ENDIF
-       ENDIF
-       IF ( nest_bound_n )  THEN   !  North border pe
-          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
-                                    nzb_u_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'u' )
-          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
-                                    nzb_v_inner, logc_v_n, logc_ratio_v_n, nzt_topo_nestbc_n, 'n', 'v' )
-          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
-                                    nzb_w_inner, logc_w_n, logc_ratio_w_n, nzt_topo_nestbc_n, 'n', 'w' )
-          CALL pmci_interp_tril_sn( e,  ec,  ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                  &
-                                    nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'e' )
-          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
-                                    nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
+
+       ENDIF
+!
+!--    North border pe
+       IF ( nest_bound_n )  THEN
+          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_n,   &
+                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
+                                    'u' )
+          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
+                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_n,   &
+                                    logc_ratio_v_n, nzt_topo_nestbc_n, 'n',    &
+                                    'v' )
+          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_n,   &
+                                    logc_ratio_w_n, nzt_topo_nestbc_n, 'n',    &
+                                    'w' )
+          CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_n,   &
+                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
+                                    'e' )
+          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
+                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_n,   &
+                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
+                                    's' )
           IF ( humidity  .OR.  passive_scalar )  THEN
-             CALL pmci_interp_tril_sn( q, qc, ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                 &
-                                       nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
+             CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
+                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
+                                       logc_u_n, logc_ratio_u_n,               &
+                                       nzt_topo_nestbc_n, 'n', 's' )
           ENDIF
+
           IF ( nesting_mode == 'one-way' )  THEN
              CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 'n', 'u' )
@@ -2518,17 +3187,26 @@
                 CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 'n', 's' )
              ENDIF
+
           ENDIF
-      ENDIF
-
-!
-!--   All PEs are top-border PEs
-       CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, 'u' )
-       CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, 'v' )
-       CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, 'w' )
-       CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 'e' )
-       CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
-       IF ( humidity  .OR.  passive_scalar )  THEN
-          CALL pmci_interp_tril_t( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
-       ENDIF
+
+       ENDIF
+
+!
+!--    All PEs are top-border PEs
+       CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
+                                r2yo, r1zo, r2zo, 'u' )
+       CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,      &
+                                r2yv, r1zo, r2zo, 'v' )
+       CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
+                                r2yo, r1zw, r2zw, 'w' )
+       CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,      &
+                                r2yo, r1zo, r2zo, 'e' )
+       CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,      &
+                                r2yo, r1zo, r2zo, 's' )
+       IF ( humidity .OR. passive_scalar )  THEN
+          CALL pmci_interp_tril_t( q, qc, ico, jco, kco, r1xo, r2xo, r1yo,     &
+                                   r2yo, r1zo, r2zo, 's' )
+       ENDIF
+
        IF ( nesting_mode == 'one-way' )  THEN
           CALL pmci_extrap_ifoutflow_t( u,  'u' )
@@ -2541,4 +3219,5 @@
           ENDIF
       ENDIF
+
    END SUBROUTINE pmci_interpolation
 
@@ -2551,10 +3230,15 @@
       IMPLICIT NONE
 
-      CALL pmci_anterp_tophat( u,  uc,  kceu, iflu, ifuu, jflo, jfuo, kflo, kfuo, nzb_u_inner, 'u' )
-      CALL pmci_anterp_tophat( v,  vc,  kceu, iflo, ifuo, jflv, jfuv, kflo, kfuo, nzb_v_inner, 'v' )
-      CALL pmci_anterp_tophat( w,  wc,  kcew, iflo, ifuo, jflo, jfuo, kflw, kfuw, nzb_w_inner, 'w' )
-      CALL pmci_anterp_tophat( pt, ptc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
+      CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
+                               kfuo, nzb_u_inner, 'u' )
+      CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
+                               kfuo, nzb_v_inner, 'v' )
+      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
+                               kfuw, nzb_w_inner, 'w' )
+      CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo,    &
+                               kfuo, nzb_s_inner, 's' )
       IF ( humidity  .OR.  passive_scalar )  THEN
-         CALL pmci_anterp_tophat( q, qc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
+         CALL pmci_anterp_tophat( q, qc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
+                                  kfuo, nzb_s_inner, 's' )
       ENDIF
 
@@ -2563,34 +3247,16 @@
 
 
-   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, &
-                                  nzt_topo_nestbc, edge, var )
-
+   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
+                                   r2z, kb, logc, logc_ratio, nzt_topo_nestbc, &
+                                   edge, var )
 !
 !--   Interpolation of ghost-node values used as the client-domain boundary
-!--   conditions. This subroutine handles the left and right boundaries. 
-!--   This subroutine is based on trilinear interpolation.
-!--   Constant dz is still assumed.
-!
-!--                                   Antti Hellsten 22.2.2015.
-!
-!--   Rewritten so that all the coefficients and client-array indices are 
-!--   precomputed in the initialization phase by pmci_init_interp_tril.
-!
-!                                     Antti Hellsten 3.3.2015.
-!
-!--   Constant dz no more assumed.
-!                                     Antti Hellsten 23.3.2015.
-!
-!--   Adapted for non-flat topography. However, the near-wall velocities 
-!--   are log-corrected only over horizontal surfaces, not yet near vertical 
-!--   walls.
-!--                                   Antti Hellsten 26.3.2015.
-!
-!--   Indexing in the principal direction (i) is changed. Now, the nest-boundary 
-!--   values are interpolated only into the first ghost-node layers on each later 
-!--   boundary. These values are then simply copied to the second ghost-node layer.
-!
-!--                                   Antti Hellsten 6.10.2015.
+!--   conditions. This subroutine handles the left and right boundaries. It is
+!--   based on trilinear interpolation. Constant dz is still assumed.
+!--   TO_DO:  constant dz is an important restriction and should be checked
+!--           somewhere in order to let users not run into this trap
       IMPLICIT NONE
+
+!--   TO_DO: wrap long lines in this and the remaining interp_tril routines
       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT)          ::  f            !:
       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)                 ::  fc           !:
@@ -2615,4 +3281,5 @@
       INTEGER(iwp) ::  i       !:
       INTEGER(iwp) ::  ib      !:
+      INTEGER(iwp) ::  ibgp    !:
       INTEGER(iwp) ::  iw      !:
       INTEGER(iwp) ::  j       !:
@@ -2759,6 +3426,14 @@
 
 !
-!--   Store the boundary values also into the second ghost node layer.
-      f(0:nzt+1,nys:nyn+1,ib) = f(0:nzt+1,nys:nyn+1,i)
+!--   Store the boundary values also into the other redundant ghost node layers
+      IF ( edge == 'l' )  THEN
+         DO ibgp = -nbgp, ib
+            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
+         ENDDO
+      ELSE IF ( edge == 'r' )  THEN
+         DO ibgp = ib, nx + nbgp
+            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
+         ENDDO
+      ENDIF
 
    END SUBROUTINE pmci_interp_tril_lr
@@ -2823,4 +3498,5 @@
       INTEGER(iwp) ::  j       !:
       INTEGER(iwp) ::  jb      !:
+      INTEGER(iwp) ::  jbgp    !:
       INTEGER(iwp) ::  k       !:
       INTEGER(iwp) ::  kcorr   !:
@@ -2960,6 +3636,14 @@
 
 !
-!--   Store the boundary values also into the second ghost node layer.
-      f(0:nzt+1,jb,nxl:nxr+1) = f(0:nzt+1,j,nxl:nxr+1)
+!--   Store the boundary values also into the other redundant ghost node layers
+      IF ( edge == 's' )  THEN
+         DO jbgp = -nbgp, jb
+            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
+         ENDDO
+      ELSE IF ( edge == 'n' )  THEN
+         DO jbgp = jb, ny + nbgp
+            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
+         ENDDO
+      ENDIF
 
    END SUBROUTINE pmci_interp_tril_sn
@@ -2974,22 +3658,4 @@
 !--   This subroutine is based on trilinear interpolation.
 !--   Constant dz is still assumed.
-!
-!--                                      Antti Hellsten 23.2.2015.
-!
-!
-!--   Rewritten so that all the coefficients and client-array indices are 
-!--   precomputed in the initialization phase by pmci_init_interp_tril.
-!
-!--                                      Antti Hellsten 3.3.2015.
-!
-!--   Constant dz no more assumed.
-!--                                      Antti Hellsten 23.3.2015.
-!
-!--   Indexing in the principal direction (k) is changed. Now, the nest-boundary 
-!--   values are interpolated only into the first ghost-node layer. Actually there is 
-!--   only one ghost-node layer in the k-direction.  
-!
-!--                                      Antti Hellsten 6.10.2015.
-!
       IMPLICIT NONE
       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
@@ -3069,839 +3735,327 @@
 
 
-   SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
-
-!
-!--   After the interpolation of ghost-node values for the client-domain boundary
-!--   conditions, this subroutine checks if there is a local outflow through the 
-!--   boundary. In that case this subroutine overwrites the interpolated values 
-!--   by values extrapolated from the domain. This subroutine handles the left and 
-!--   right boundaries.
-!--   However, this operation is only needed in case of one-way coupling.
-!
-!--                                       Antti Hellsten 9.3.2015.
-!
-!--   Indexing in the principal direction (i) is changed. Now, the nest-boundary 
-!--   values are interpolated only into the first ghost-node layers on each later 
-!--   boundary. These values are then simply copied to the second ghost-node layer.
-!
-!--                                       Antti Hellsten 6.10.2015.
-!
-      IMPLICIT NONE
-      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
-
-      CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
-      CHARACTER(LEN=1),INTENT(IN) ::  var    !:
-
-      INTEGER(iwp) ::  i     !:
-      INTEGER(iwp) ::  ib    !:
-      INTEGER(iwp) ::  ied   !:
-      INTEGER(iwp) ::  j     !:
-      INTEGER(iwp) ::  k     !:
+    SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
+!
+!--    After the interpolation of ghost-node values for the client-domain
+!--    boundary conditions, this subroutine checks if there is a local outflow
+!--    through the boundary. In that case this subroutine overwrites the
+!--    interpolated values by values extrapolated from the domain. This
+!--    subroutine handles the left and right boundaries. However, this operation
+!--    is only needed in case of one-way coupling.
+       IMPLICIT NONE
+
+       CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
+       CHARACTER(LEN=1),INTENT(IN) ::  var    !:
+
+       INTEGER(iwp) ::  i     !:
+       INTEGER(iwp) ::  ib    !:
+       INTEGER(iwp) ::  ibgp  !:
+       INTEGER(iwp) ::  ied   !:
+       INTEGER(iwp) ::  j     !:
+       INTEGER(iwp) ::  k     !:
       
-      REAL(wp) ::  outnor    !:
-      REAL(wp) ::  vdotnor   !:
-
-!
-!--   Check which edge is to be handled: left or right. 
-      IF ( edge == 'l' )  THEN
-         IF ( var == 'u' )  THEN
-            i   = nxl
-            ib  = nxl - 1
-            ied = nxl + 1
-         ELSE
-            i   = nxl - 1
-            ib  = nxl - 2
-            ied = nxl
-         ENDIF
-         outnor = -1.0_wp
-      ELSEIF ( edge == 'r' )  THEN
-         i      = nxr + 1
-         ib     = nxr + 2
-         ied    = nxr
-         outnor = 1.0_wp
-      ENDIF
-
-      DO  j = nys, nyn + 1
-         DO  k = kb(j,i), nzt +1
-            vdotnor = outnor * u(k,j,ied)
-            IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
-               f(k,j,i) = f(k,j,ied)
-            ENDIF
-         ENDDO
-         IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
-            f(kb(j,i),j,i) = 0.0_wp
-         ENDIF
-      ENDDO
-
-!
-!--   Store the updated boundary values also into the second ghost node layer.
-      f(0:nzt,nys:nyn+1,ib) = f(0:nzt,nys:nyn+1,i)
-
-   END SUBROUTINE pmci_extrap_ifoutflow_lr
-
-
-
-   SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
-!
-!--   After  the interpolation of ghost-node values for the client-domain boundary
-!--   conditions, this subroutine checks if there is a local outflow through the 
-!--   boundary. In that case this subroutine overwrites the interpolated values 
-!--   by values extrapolated from the domain. This subroutine handles the south and 
-!--   north boundaries. 
-!
-!--                                       Antti Hellsten 9.3.2015.
-!
-!--   Indexing in the principal direction (j) is changed. Now, the nest-boundary 
-!--   values are interpolated only into the first ghost-node layers on each later 
-!--   boundary. These values are then simply copied to the second ghost-node layer.
-!
-!--                                       Antti Hellsten 6.10.2015.
-
-      IMPLICIT NONE
-      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
-      CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
-      CHARACTER(LEN=1), INTENT(IN) ::  var    !:
+       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
+
+       REAL(wp) ::  outnor    !:
+       REAL(wp) ::  vdotnor   !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
+
+!
+!--    Check which edge is to be handled: left or right
+       IF ( edge == 'l' )  THEN
+          IF ( var == 'u' )  THEN
+             i   = nxl
+             ib  = nxl - 1
+             ied = nxl + 1
+          ELSE
+             i   = nxl - 1
+             ib  = nxl - 2
+             ied = nxl
+          ENDIF
+          outnor = -1.0_wp
+       ELSEIF ( edge == 'r' )  THEN
+          i      = nxr + 1
+          ib     = nxr + 2
+          ied    = nxr
+          outnor = 1.0_wp
+       ENDIF
+
+       DO  j = nys, nyn + 1
+          DO  k = kb(j,i), nzt +1
+             vdotnor = outnor * u(k,j,ied)
+!
+!--          Local outflow
+             IF ( vdotnor > 0.0_wp )  THEN
+                f(k,j,i) = f(k,j,ied)
+             ENDIF
+          ENDDO
+          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
+             f(kb(j,i),j,i) = 0.0_wp
+          ENDIF
+       ENDDO
+
+!
+!--    Store the boundary values also into the redundant ghost node layers.
+       IF ( edge == 'l' )  THEN
+          DO ibgp = -nbgp, ib
+             f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
+          ENDDO
+       ELSEIF ( edge == 'r' )  THEN
+          DO ibgp = ib, nx + nbgp
+             f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
+          ENDDO
+       ENDIF
+
+    END SUBROUTINE pmci_extrap_ifoutflow_lr
+
+
+
+    SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
+!
+!--    After  the interpolation of ghost-node values for the client-domain
+!--    boundary conditions, this subroutine checks if there is a local outflow
+!--    through the boundary. In that case this subroutine overwrites the
+!--    interpolated values by values extrapolated from the domain. This
+!--    subroutine handles the south and north boundaries.
+       IMPLICIT NONE
+
+       CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
+       CHARACTER(LEN=1), INTENT(IN) ::  var    !:
       
-      INTEGER(iwp) ::  i         !:
-      INTEGER(iwp) ::  j         !:
-      INTEGER(iwp) ::  jb        !:
-      INTEGER(iwp) ::  jed       !:
-      INTEGER(iwp) ::  k         !:
-      REAL(wp)     ::  outnor    !:
-      REAL(wp)     ::  vdotnor   !:
-
-!
-!--   Check which edge is to be handled: left or right. 
-      IF ( edge == 's' )  THEN
-         IF ( var == 'v' )  THEN
-            j   = nys
-            jb  = nys - 1
-            jed = nys + 1
-         ELSE
-            j   = nys - 1
-            jb  = nys - 2
-            jed = nys
-         ENDIF
-         outnor = -1.0_wp
-      ELSEIF ( edge == 'n' )  THEN
-         j      = nyn + 1
-         jb     = nyn + 2
-         jed    = nyn
-         outnor = 1.0_wp
-      ENDIF
-
-      DO  i = nxl, nxr + 1
-         DO  k = kb(j,i), nzt + 1
-            vdotnor = outnor * v(k,jed,i)
-            IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
-               f(k,j,i) = f(k,jed,i)                 
-            ENDIF
-         ENDDO
-         IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
-            f(kb(j,i),j,i) = 0.0_wp
-         ENDIF
-      ENDDO
-
-!
-!--   Store the updated boundary values also into the second ghost node layer.
-      f(0:nzt,jb,nxl:nxr+1) = f(0:nzt,j,nxl:nxr+1)
-
-   END SUBROUTINE pmci_extrap_ifoutflow_sn
+       INTEGER(iwp) ::  i         !:
+       INTEGER(iwp) ::  j         !:
+       INTEGER(iwp) ::  jb        !:
+       INTEGER(iwp) ::  jbgp      !:
+       INTEGER(iwp) ::  jed       !:
+       INTEGER(iwp) ::  k         !:
+
+       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
+
+       REAL(wp)     ::  outnor    !:
+       REAL(wp)     ::  vdotnor   !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
+
+!
+!--    Check which edge is to be handled: left or right
+       IF ( edge == 's' )  THEN
+          IF ( var == 'v' )  THEN
+             j   = nys
+             jb  = nys - 1
+             jed = nys + 1
+          ELSE
+             j   = nys - 1
+             jb  = nys - 2
+             jed = nys
+          ENDIF
+          outnor = -1.0_wp
+       ELSEIF ( edge == 'n' )  THEN
+          j      = nyn + 1
+          jb     = nyn + 2
+          jed    = nyn
+          outnor = 1.0_wp
+       ENDIF
+
+       DO  i = nxl, nxr + 1
+          DO  k = kb(j,i), nzt + 1
+             vdotnor = outnor * v(k,jed,i)
+!
+!--          Local outflow
+             IF ( vdotnor > 0.0_wp )  THEN
+                f(k,j,i) = f(k,jed,i)
+             ENDIF
+          ENDDO
+          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
+             f(kb(j,i),j,i) = 0.0_wp
+          ENDIF
+       ENDDO
+
+!
+!--    Store the boundary values also into the redundant ghost node layers.
+       IF ( edge == 's' )  THEN
+          DO jbgp = -nbgp, jb
+             f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
+          ENDDO
+       ELSEIF ( edge == 'n' )  THEN
+          DO jbgp = jb, ny + nbgp
+             f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
+          ENDDO
+       ENDIF
+
+    END SUBROUTINE pmci_extrap_ifoutflow_sn
 
  
 
-   SUBROUTINE pmci_extrap_ifoutflow_t( f, var )
-
-!
-!--   Interpolation of ghost-node values used as the client-domain boundary
-!--   conditions. This subroutine handles the top boundary. 
-!--   This subroutine is based on trilinear interpolation.
-!     
-!--                                       Antti Hellsten 23.2.2015.
-!     
-!     
-!--   Rewritten so that all the coefficients and client-array indices are 
-!--   precomputed in the initialization phase by init_interp_tril.
-!     
-!--                                       Antti Hellsten 3.3.2015. 
-!     
-!--   Indexing in the principal direction (k) is changed. Now, the nest-boundary 
-!--   values are extrapolated only into the first ghost-node layer. Actually there is 
-!--   only one ghost-node layer in the k-direction.  
-!     
-!--                                       Antti Hellsten 6.10.2015.
-      IMPLICIT NONE
-      REAL(wp), DIMENSION(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp), INTENT(INOUT) ::  f   !:
-      CHARACTER(LEN=1), INTENT(IN) ::  var   !:
+    SUBROUTINE pmci_extrap_ifoutflow_t( f, var )
+!
+!--    Interpolation of ghost-node values used as the client-domain boundary
+!--    conditions. This subroutine handles the top boundary. It is based on
+!--    trilinear interpolation.
+       IMPLICIT NONE
+
+       CHARACTER(LEN=1), INTENT(IN) ::  var   !:
       
-      INTEGER(iwp) ::  i     !:
-      INTEGER(iwp) ::  j     !:
-      INTEGER(iwp) ::  k     !:
-      INTEGER(iwp) ::  ked   !:
-      REAL(wp) ::  vdotnor   !:
+       INTEGER(iwp) ::  i     !:
+       INTEGER(iwp) ::  j     !:
+       INTEGER(iwp) ::  k     !:
+       INTEGER(iwp) ::  ked   !:
+
+       REAL(wp) ::  vdotnor   !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp),     &
+                 INTENT(INOUT) ::  f   !:
      
 
-      IF ( var == 'w' )  THEN
-         k    = nzt
-         ked  = nzt - 1
-      ELSE
-         k    = nzt + 1  
-         ked  = nzt
-      ENDIF
-
-      DO  i = nxl, nxr
-         DO  j = nys, nyn
-            vdotnor = w(ked,j,i)
-            IF ( vdotnor > 0.0_wp )  THEN   !:  Local outflow.              
-               f(k,j,i) = f(ked,j,i)
-            ENDIF
-         ENDDO
-      ENDDO
-
-!
-!--   Just fill up the second ghost-node layer for w.
-      IF ( var == 'w' )  THEN
-         f(nzt+1,:,:) = f(nzt,:,:)
-      ENDIF
-
-   END SUBROUTINE pmci_extrap_ifoutflow_t
-
-
-
-   SUBROUTINE pmci_anterp_tophat( f, fc, kce, ifl, ifu, jfl, jfu, kfl, kfu, kb, var )
-!
-!--   Anterpolation of internal-node values to be used as the server-domain 
-!--   values. This subroutine is based on the first-order numerical 
-!--   integration of the fine-grid values contained within the coarse-grid
-!--   cell. 
-!     
-!--                                        Antti Hellsten 16.9.2015.
-!
-      IMPLICIT NONE
-      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN)  ::  f     !:
-      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)   ::  fc    !:
-      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifl   !:
-      INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifu   !:
-      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfl   !:
-      INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfu   !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)        ::  kb    !: may be unnecessary
-      INTEGER(iwp), INTENT(IN)                                        ::  kce   !:
-      INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfl   !:
-      INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfu   !:
-      CHARACTER(LEN=1), INTENT(IN)                                    ::  var   !:                            
-
-      INTEGER(iwp) ::  i         !:
-      INTEGER(iwp) ::  icb       !:
-      INTEGER(iwp) ::  ice       !:
-      INTEGER(iwp) ::  ifc       !:
-      INTEGER(iwp) ::  ijfc      !:
-      INTEGER(iwp) ::  j         !:
-      INTEGER(iwp) ::  jcb       !:
-      INTEGER(iwp) ::  jce       !:
-      INTEGER(iwp) ::  k         !:
-      INTEGER(iwp) ::  kcb       !:
-      INTEGER(iwp) ::  l         !:
-      INTEGER(iwp) ::  m         !:
-      INTEGER(iwp) ::  n         !:
-      INTEGER(iwp) ::  nfc       !:
-      REAL(wp)     ::  cellsum   !:
-      REAL(wp)     ::  f1f       !:
-      REAL(wp)     ::  fra       !:
+       IF ( var == 'w' )  THEN
+          k    = nzt
+          ked  = nzt - 1
+       ELSE
+          k    = nzt + 1
+          ked  = nzt
+       ENDIF
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             vdotnor = w(ked,j,i)
+!
+!--          Local outflow
+             IF ( vdotnor > 0.0_wp )  THEN
+                f(k,j,i) = f(ked,j,i)
+             ENDIF
+          ENDDO
+       ENDDO
+
+!
+!--    Just fill up the second ghost-node layer for w
+       IF ( var == 'w' )  THEN
+          f(nzt+1,:,:) = f(nzt,:,:)
+       ENDIF
+
+    END SUBROUTINE pmci_extrap_ifoutflow_t
+
+
+
+    SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu,   &
+                                   kb, var )
+!
+!--    Anterpolation of internal-node values to be used as the server-domain
+!--    values. This subroutine is based on the first-order numerical
+!--    integration of the fine-grid values contained within the coarse-grid
+!--    cell.
+       IMPLICIT NONE
+
+       CHARACTER(LEN=1), INTENT(IN) ::  var   !:
+
+       INTEGER(iwp) ::  i         !: Fine-grid index
+       INTEGER(iwp) ::  ii        !: Coarse-grid index
+       INTEGER(iwp) ::  iclp      !:
+       INTEGER(iwp) ::  icrm      !:
+       INTEGER(iwp) ::  ifc       !:
+       INTEGER(iwp) ::  ijfc      !:
+       INTEGER(iwp) ::  j         !: Fine-grid index
+       INTEGER(iwp) ::  jj        !: Coarse-grid index
+       INTEGER(iwp) ::  jcnm      !:
+       INTEGER(iwp) ::  jcsp      !:
+       INTEGER(iwp) ::  k         !: Fine-grid index
+       INTEGER(iwp) ::  kk        !: Coarse-grid index
+       INTEGER(iwp) ::  kcb       !:
+       INTEGER(iwp) ::  nfc       !:
+
+       INTEGER(iwp), INTENT(IN) ::  kct   !:
+
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)             ::  ifl   !:
+       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)             ::  ifu   !:
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)             ::  jfl   !:
+       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)             ::  jfu   !:
+!--    TO_DO: is the next line really unnecessary?
+       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !: may be unnecessary
+       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)               ::  kfl   !:
+       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)               ::  kfu   !:
+
+
+       REAL(wp) ::  cellsum   !:
+       REAL(wp) ::  f1f       !:
+       REAL(wp) ::  fra       !:
+
+       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  f   !:
+       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)  ::  fc  !:
   
 
-      icb = icl 
-      ice = icr 
-      jcb = jcs 
-      jce = jcn 
-
-!
-!--   Define the index bounds icb, ice, jcb and jce. 
-!--   Note that kcb is simply zero and kce enters here as a parameter and it is 
-!--   determined in init_anterp_tophat  
-      IF ( nest_bound_l )  THEN
-         IF ( var == 'u' )  THEN
-            icb = icl + nhll + 1
-         ELSE
-            icb = icl + nhll
-         ENDIF
-      ENDIF
-      IF ( nest_bound_r )  THEN
-         ice = icr - nhlr
-      ENDIF
-
-      IF ( nest_bound_s )  THEN
-         IF ( var == 'v' )  THEN
-            jcb = jcs + nhls + 1
-         ELSE
-            jcb = jcs + nhls
-         ENDIF
-      ENDIF
-      IF ( nest_bound_n )  THEN
-         jce = jcn - nhln
-      ENDIF
-      kcb = 0 
-
-!
-!--   Note that l,m, and n are coarse-grid indices and i,j, and k are fine-grid indices.
-      DO  l = icb, ice
-         ifc = ifu(l) - ifl(l) + 1
-         DO  m = jcb, jce
-            ijfc = ifc * ( jfu(m) - jfl(m) +1 )
-
-!
-!--         How to deal with the lower bound of k in case of non-flat topography?
-            !kcb = MIN( kb(jfl(m),ifl(l)), kb(jfu(m),ifl(l)), kb(jfl(m),ifu(l)), kb(jfu(m),ifu(l)) ) ! Something wrong with this.
-            DO  n = kcb, kce
-               nfc =  ijfc * ( kfu(n) - kfl(n) + 1 )
-               cellsum = 0.0
-               DO  i = ifl(l), ifu(l)
-                  DO  j = jfl(m), jfu(m)
-                     DO  k = kfl(n), kfu(n)
-                        cellsum = cellsum + f(k,j,i)
-                     ENDDO
-                  ENDDO
-               ENDDO
-
-!
-!--            Spatial under-relaxation.
-               fra  = frax(l) * fray(m) * fraz(n)
-!--            TO_DO: why not KIND=wp ?
-               fc(n,m,l) = ( 1.0_wp - fra ) * fc(n,m,l) + fra * cellsum / REAL( nfc, KIND=KIND(cellsum) )  
-            ENDDO
-         ENDDO
-      ENDDO
-
-   END SUBROUTINE pmci_anterp_tophat
+!
+!--    Initialize the index bounds for anterpolation 
+       iclp = icl 
+       icrm = icr 
+       jcsp = jcs 
+       jcnm = jcn 
+
+!
+!--    Define the index bounds iclp, icrm, jcsp and jcnm.
+!--    Note that kcb is simply zero and kct enters here as a parameter and it is
+!--    determined in pmci_init_anterp_tophat
+       IF ( nest_bound_l )  THEN
+          IF ( var == 'u' )  THEN
+             iclp = icl + nhll + 1
+          ELSE
+             iclp = icl + nhll
+          ENDIF
+       ENDIF
+       IF ( nest_bound_r )  THEN
+          icrm = icr - nhlr
+       ENDIF
+
+       IF ( nest_bound_s )  THEN
+          IF ( var == 'v' )  THEN
+             jcsp = jcs + nhls + 1
+          ELSE
+             jcsp = jcs + nhls
+          ENDIF
+       ENDIF
+       IF ( nest_bound_n )  THEN
+          jcnm = jcn - nhln
+       ENDIF
+       kcb = 0
+
+!
+!--    Note that l,m, and n are coarse-grid indices and i,j, and k are fine-grid
+!--    indices.
+       DO  ii = iclp, icrm
+          ifc = ifu(ii) - ifl(ii) + 1
+          DO  jj = jcsp, jcnm
+             ijfc = ifc * ( jfu(jj) - jfl(jj) + 1 )
+!
+!--          For simplicity anterpolate within buildings too
+             DO  kk = kcb, kct
+                nfc =  ijfc * ( kfu(kk) - kfl(kk) + 1 )
+                cellsum = 0.0_wp
+                DO  i = ifl(ii), ifu(ii)
+                   DO  j = jfl(jj), jfu(jj)
+                      DO  k = kfl(kk), kfu(kk)
+                         cellsum = cellsum + f(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+!
+!--             Spatial under-relaxation.
+                fra  = frax(ii) * fray(jj) * fraz(kk)
+!
+!--             Block out the fine-grid corner patches from the anterpolation
+                IF ( ( ifl(ii) < nxl ) .OR. ( ifu(ii) > nxr ) )  THEN
+                   IF ( ( jfl(jj) < nys ) .OR. ( jfu(jj) > nyn ) )  THEN
+                      fra = 0.0_wp
+                   ENDIF
+                ENDIF
+
+                fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +                 &
+                               fra * cellsum / REAL( nfc, KIND = wp )
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE pmci_anterp_tophat
 
 #endif
  END SUBROUTINE pmci_client_datatrans
 
-
-
- SUBROUTINE pmci_update_new
-
-#if defined( __parallel )
-!
-!-- Copy the interpolated/anterpolated boundary values to the _p 
-!-- arrays, too, to make sure the interpolated/anterpolated boundary
-!-- values are carried over from one RK inner step to another.
-!-- So far works only with the cpp-switch __nopointer.
-!
-!--                 Antti Hellsten 8.3.2015 
-!
-
-!-- Just debugging
-    w(nzt+1,:,:) = w(nzt,:,:)
-
-    u_p  = u
-    v_p  = v
-    w_p  = w
-    e_p  = e
-    pt_p = pt
-    IF ( humidity  .OR.  passive_scalar )  THEN
-       q_p = q
-    ENDIF
-
-!
-!-- TO_DO: Find out later if nesting would work without __nopointer.
-#endif
-
- END SUBROUTINE pmci_update_new
-
-
-
- SUBROUTINE pmci_set_array_pointer( name, client_id, nz_cl )
-
-    IMPLICIT NONE
-
-    INTEGER, INTENT(IN)                  ::  client_id   !:
-    INTEGER, INTENT(IN)                  ::  nz_cl       !:
-    CHARACTER(LEN=*), INTENT(IN)         ::  name        !:
-    
-#if defined( __parallel )
-    REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d        !:
-    REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d_sec    !:
-    REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d        !:
-    REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d_sec    !:
-    INTEGER(iwp)                         ::  ierr        !:
-    INTEGER(iwp)                         ::  istat       !:
-    
-
-    NULLIFY( p_3d )
-    NULLIFY( p_2d )
-
-!
-!-- List of array names, which can be coupled
-!-- In case of 3D please change also the second array for the pointer version
-    IF ( TRIM(name) == "u" )     p_3d => u
-    IF ( TRIM(name) == "v" )     p_3d => v
-    IF ( TRIM(name) == "w" )     p_3d => w
-    IF ( TRIM(name) == "e" )     p_3d => e
-    IF ( TRIM(name) == "pt" )    p_3d => pt
-    IF ( TRIM(name) == "q" )     p_3d => q
-!
-!-- This 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 defined( __nopointer )
-    IF ( ASSOCIATED( p_3d ) )  THEN
-       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz )
-    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
-       CALL pmc_s_set_dataarray( client_id, p_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'
-          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
-       ELSE
-!
-!--       Avoid others to continue
-          CALL MPI_BARRIER( comm2d, ierr )
-       ENDIF
-    ENDIF
-#else
-!-- TO_DO: Why aren't the other pointers (p_3d) not set to u_1, v_1, etc.??
-    IF ( TRIM(name) == "u" )     p_3d_sec => u_2
-    IF ( TRIM(name) == "v" )     p_3d_sec => v_2
-    IF ( TRIM(name) == "w" )     p_3d_sec => w_2
-    IF ( TRIM(name) == "e" )     p_3d_sec => e_2
-    IF ( TRIM(name) == "pt" )    p_3d_sec => pt_2
-    IF ( TRIM(name) == "q" )     p_3d_sec => q_2
-
-    IF ( ASSOCIATED( p_3d ) )  THEN
-       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz, &
-                                 array_2 = p_3d_sec )
-    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
-       CALL pmc_s_set_dataarray( client_id, p_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'
-          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
-       ELSE
-!
-!--       Avoid others to continue
-          CALL MPI_BARRIER( comm2d, ierr )
-       ENDIF
-
-   ENDIF
-#endif
-
-#endif
- END SUBROUTINE pmci_set_array_pointer
-
-
-
- SUBROUTINE pmci_create_client_arrays( name, is, ie, js, je, nzc  )
-
-    IMPLICIT NONE
-
-    INTEGER(iwp), INTENT(IN)              ::  ie      !:
-    INTEGER(iwp), INTENT(IN)              ::  is      !:
-    INTEGER(iwp), INTENT(IN)              ::  je      !:
-    INTEGER(iwp), INTENT(IN)              ::  js      !:
-    INTEGER(iwp), INTENT(IN)              ::  nzc     !:  Note that nzc is cg%nz
-    CHARACTER(LEN=*), INTENT(IN)          ::  name    !:
-     
-#if defined( __parallel )
-    REAL(wp), POINTER,DIMENSION(:,:)      ::  p_2d    !:
-    REAL(wp), POINTER,DIMENSION(:,:,:)    ::  p_3d    !:
-    INTEGER(iwp)                          ::  ierr    !:
-    INTEGER(iwp)                          ::  istat   !:
-
-
-    NULLIFY( p_3d )
-    NULLIFY( p_2d )
-
-!
-!-- List of array names, which can be coupled.
-!-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1.
-    IF ( TRIM( name ) == "u" )  THEN
-       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
-       p_3d => uc
-    ELSEIF ( TRIM( name ) == "v" )  THEN
-       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1, js:je, is:ie) )
-       p_3d => vc
-    ELSEIF ( TRIM( name ) == "w" )  THEN
-       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1, js:je, is:ie) )
-       p_3d => wc
-    ELSEIF ( TRIM( name ) == "e" )  THEN
-       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1, js:je, is:ie) )
-       p_3d => ec
-    ELSEIF ( TRIM( name ) == "pt")  THEN
-       IF ( .NOT. ALLOCATED( ptc ) ) ALLOCATE( ptc(0:nzc+1, js:je, is:ie) )
-       p_3d => ptc
-    ELSEIF ( TRIM( name ) == "q")  THEN
-       IF ( .NOT. ALLOCATED( qc ) ) ALLOCATE( qc(0:nzc+1, js:je, is:ie) )
-       p_3d => qc
-    !ELSEIF (trim(name) == "z0") then
-       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
-       !p_2d => z0c
-    ENDIF
-
-    IF ( ASSOCIATED( p_3d ) )  THEN
-       CALL pmc_c_set_dataarray( p_3d )
-    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
-       CALL pmc_c_set_dataarray( p_2d )
-    ELSE
-!
-!--    Give only one message for the first client domain
-       IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
-
-          message_string = 'pointer for array "' // TRIM( name ) //            &
-                           '" can''t be associated'
-          CALL message( 'pmci_create_client_arrays', 'PA0170', 3, 2, 0, 6, 0 )
-       ELSE
-!
-!--       Avoid others to continue
-          CALL MPI_BARRIER( comm2d, ierr )
-       ENDIF
-    ENDIF
-
-#endif
- END SUBROUTINE pmci_create_client_arrays
-
-
-
- SUBROUTINE pmci_server_initialize
-
-#if defined( __parallel )
-    IMPLICIT NONE
-
-    INTEGER(iwp)   ::  client_id   !:
-    INTEGER(iwp)   ::  m           !:
-    REAL(wp)       ::  waittime    !:
-
-
-    DO  m = 1, SIZE( pmc_server_for_client ) - 1
-       client_id = pmc_server_for_client(m)
-       CALL pmc_s_fillbuffer( client_id, waittime=waittime )
-    ENDDO
-
-#endif
- END SUBROUTINE pmci_server_initialize
-
-
-
- SUBROUTINE pmci_client_initialize
-
-#if defined( __parallel )
-    IMPLICIT NONE
-
-    INTEGER(iwp)   ::  i          !:
-    INTEGER(iwp)   ::  icl        !:
-    INTEGER(iwp)   ::  icr        !:
-    INTEGER(iwp)   ::  j          !:
-    INTEGER(iwp)   ::  jcn        !:
-    INTEGER(iwp)   ::  jcs        !:
-    REAL(wp)       ::  waittime   !:
-
-
-    IF ( cpl_id > 1 )  THEN   !  Root id is never a client      
-
-!
-!--    Client domain boundaries in the server indice space.
-       icl = coarse_bound(1)
-       icr = coarse_bound(2)
-       jcs = coarse_bound(3)
-       jcn = coarse_bound(4)
-
-!
-!--    Get data from server      
-       CALL pmc_c_getbuffer( waittime = waittime )
-
-!
-!--    The interpolation.
-       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
-       CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
-       CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
-       CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
-       CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
-       IF ( humidity  .OR.  passive_scalar )  THEN
-          CALL pmci_interp_tril_all ( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
-       ENDIF
-
-       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
-                u(nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
-                v(nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
-                w(nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
-                e(nzb:nzb_s_inner(j,i),j,i)   = 0.0_wp
-                u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
-                v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
-                w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
-                e_p(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
-             ENDDO
-          ENDDO
-       ENDIF
-    ENDIF
-    
-
- CONTAINS
-
-
-   SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, var )
-
-!
-!-- Interpolation of the internal values for the client-domain initialization.
-!-- This subroutine is based on trilinear interpolation.
-!-- Coding based on interp_tril_lr/sn/t
-!
-!--                                     Antti Hellsten 20.10.2015.
-      IMPLICIT NONE
-      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
-      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)        ::  fc    !:
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r1x   !:
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r2x   !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r1y   !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r2y   !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r1z   !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r2z   !:
-      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                    ::  ic    !:
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                    ::  jc    !:
-      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                    ::  kc    !: 
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb    !:
-      CHARACTER(LEN=1), INTENT(IN) :: var  !:
-      
-      INTEGER(iwp) ::  i      !:
-      INTEGER(iwp) ::  ib     !:
-      INTEGER(iwp) ::  ie     !:
-      INTEGER(iwp) ::  j      !:
-      INTEGER(iwp) ::  jb     !:
-      INTEGER(iwp) ::  je     !:
-      INTEGER(iwp) ::  k      !:
-      INTEGER(iwp) ::  k1     !:
-      INTEGER(iwp) ::  kbc    !:
-      INTEGER(iwp) ::  l      !:
-      INTEGER(iwp) ::  m      !:
-      INTEGER(iwp) ::  n      !:
-      REAL(wp) ::  fk         !:
-      REAL(wp) ::  fkj        !:
-      REAL(wp) ::  fkjp       !:
-      REAL(wp) ::  fkp        !:
-      REAL(wp) ::  fkpj       !:
-      REAL(wp) ::  fkpjp      !:
-      REAL(wp) ::  logratio   !:
-      REAL(wp) ::  logzuc1    !:
-      REAL(wp) ::  zuc1       !:
-      
-      
-      ib = nxl
-      ie = nxr
-      jb = nys   
-      je = nyn
-      IF ( nest_bound_l )  THEN
-         ib = nxl - 1 
-         IF ( var == 'u' )  THEN   !  For u, nxl is a ghost node, but not for the other variables.
-            ib = nxl
-         ENDIF
-      ENDIF
-      IF ( nest_bound_s )  THEN
-         jb = nys - 1 
-         IF ( var == 'v' )  THEN   !  For v, nys is a ghost node, but not for the other variables.
-            jb = nys
-         ENDIF
-      ENDIF
-      IF ( nest_bound_r )  THEN
-         ie = nxr + 1
-      ENDIF
-      IF ( nest_bound_n )  THEN
-         je = nyn + 1
-      ENDIF
-
-!
-!--   Trilinear interpolation.
-      DO  i = ib, ie
-         DO  j = jb, je
-            DO  k = kb(j,i), nzt + 1 
-               l = ic(i)
-               m = jc(j)
-               n = kc(k)
-               fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
-               fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
-               fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
-               fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
-               fk       = r1y(j) * fkj  + r2y(j) * fkjp
-               fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
-               f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
-            ENDDO
-         ENDDO
-      ENDDO
-
-!
-!--   Correct the interpolated values of u and v in near-wall nodes, i.e. in 
-!--   the nodes below the coarse-grid nodes with k=1. The corrction is only made 
-!--   over horizontal wall surfaces in this phase. For the nest boundary conditions,
-!--   a corresponding corrections is made for all vertical walls, too. 
-      IF ( var == 'u' .OR. var == 'v' )  THEN
-         DO  i = ib, nxr
-            DO  j = jb, nyn              
-               kbc = 1
-               DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )   !  kbc is the first coarse-grid point above the surface.
-                  kbc = kbc + 1
-               ENDDO
-               zuc1    = cg%zu(kbc)
-               k1      = kb(j,i) + 1  
-               DO  WHILE ( zu(k1) < zuc1 )
-                  k1 = k1 + 1
-               ENDDO
-               logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
-               
-               k = kb(j,i) + 1
-               DO  WHILE ( zu(k) < zuc1 ) 
-                  logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) / logzuc1 
-                  f(k,j,i) = logratio * f(k1,j,i)
-                  k  = k + 1
-               ENDDO
-               f(kb(j,i),j,i) = 0.0_wp
-            ENDDO
-         ENDDO
-      ELSEIF ( var == 'w' )  THEN
-         DO  i = ib, nxr
-            DO  j = jb, nyn
-               f(kb(j,i),j,i) = 0.0_wp
-            ENDDO
-         ENDDO
-      ENDIF
-      
-   END SUBROUTINE pmci_interp_tril_all
-
-#endif
- END SUBROUTINE pmci_client_initialize
-
-
-
- SUBROUTINE pmci_ensure_nest_mass_conservation
-
-#if defined( __parallel )
-!
-!-- Adjust the volume-flow rate through the top boundary 
-!-- so that the net volume flow through all boundaries 
-!-- of the current nest domain becomes zero.
-    IMPLICIT NONE
-
-    INTEGER(iwp)  ::  i                          !:
-    INTEGER(iwp)  ::  ierr                       !:
-    INTEGER(iwp)  ::  j                          !:
-    INTEGER(iwp)  ::  k                          !:
-    REAL(wp) ::  dxdy                            !:
-    REAL(wp) ::  innor                           !:
-    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !:
-    REAL(wp) ::  w_lt                            !:
-
-!
-!-- Sum up the volume flow through the left/right boundaries.
-    volume_flow(1)   = 0.0_wp
-    volume_flow_l(1) = 0.0_wp
-
-    IF ( nest_bound_l )  THEN
-       i = 0 
-       innor = dy
-       DO   j = nys, nyn
-          DO   k = nzb_u_inner(j,i) + 1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)            
-          ENDDO
-       ENDDO
-    ENDIF
-
-    IF ( nest_bound_r )  THEN     
-       i = nx + 1
-       innor = -dy
-       DO   j = nys, nyn
-          DO   k = nzb_u_inner(j,i) + 1, nzt
-             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)            
-          ENDDO
-       ENDDO
-    ENDIF
-      
-#if defined( __parallel )   
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
-                        MPI_SUM, comm2d, ierr )
-#else
-    volume_flow(1) = volume_flow_l(1)
-#endif
-
-!
-!-- Sum up the volume flow through the south/north boundaries.
-    volume_flow(2)   = 0.0_wp
-    volume_flow_l(2) = 0.0_wp
-
-    IF ( nest_bound_s )  THEN
-       j = 0
-       innor = dx
-       DO   i = nxl, nxr
-          DO   k = nzb_v_inner(j,i) + 1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
-          ENDDO
-       ENDDO
-    ENDIF
-
-    IF ( nest_bound_n )  THEN
-       j = ny + 1
-       innor = -dx
-       DO   i = nxl, nxr
-          DO   k = nzb_v_inner(j,i)+1, nzt
-             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
-          ENDDO
-       ENDDO
-    ENDIF
-
-#if defined( __parallel )   
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
-                        MPI_SUM, comm2d, ierr )   
-#else
-    volume_flow(2) = volume_flow_l(2)
-#endif
-      
-!
-!-- Sum up the volume flow through the top boundary.
-    volume_flow(3)   = 0.0_wp
-    volume_flow_l(3) = 0.0_wp
-    dxdy = dx * dy
-    k = nzt 
-    DO   i = nxl, nxr
-       DO   j = nys, nyn
-          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
-       ENDDO
-    ENDDO
-       
-#if defined( __parallel )   
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL, &
-                        MPI_SUM, comm2d, ierr )   
-#else
-    volume_flow(3) = volume_flow_l(3)
-#endif
-
-!
-!-- Correct the top-boundary value of w.
-    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
-    DO   i = nxl, nxr
-       DO   j = nys, nyn
-          DO  k = nzt, nzt + 1
-             w(k,j,i) = w(k,j,i) + w_lt
-          ENDDO
-       ENDDO
-    ENDDO
-
-#endif
- END SUBROUTINE pmci_ensure_nest_mass_conservation
-
 END MODULE pmc_interface
Index: palm/trunk/SOURCE/pmc_server.f90
===================================================================
--- palm/trunk/SOURCE/pmc_server.f90	(revision 1789)
+++ palm/trunk/SOURCE/pmc_server.f90	(revision 1791)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Debug write-statements commented out
 !
 ! Former revisions:
@@ -136,5 +136,5 @@
 
       do i=1,size(PMC_Server_for_Client)-1
-         if(m_model_comm == 0) write(0,*) 'PMC_Server: Initialize client Id',PMC_Server_for_Client(i)
+!         if(m_model_comm == 0) write(0,*) 'PMC_Server: Initialize client Id',PMC_Server_for_Client(i)
 
          ClientId = PMC_Server_for_Client(i)
@@ -152,5 +152,5 @@
          CALL MPI_Comm_rank (Clients(ClientId)%intra_comm, Clients(ClientId)%intra_rank, istat);
 
-         write(9,*) 'ClientId ',i,ClientId,m_world_rank, Clients(ClientId)%inter_npes
+!         write(9,*) 'ClientId ',i,ClientId,m_world_rank, Clients(ClientId)%inter_npes
 
          ALLOCATE (Clients(ClientId)%PEs(Clients(ClientId)%inter_npes))
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1789)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1791)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! call of pmci_update_new removed
 !
 ! Former revisions:
@@ -257,5 +257,5 @@
                pmci_ensure_nest_mass_conservation, pmci_client_datatrans,      &
                pmci_client_synchronize, pmci_server_datatrans,                 &
-               pmci_server_synchronize, pmci_update_new, server_to_client
+               pmci_server_synchronize, server_to_client
 
     USE production_e_mod,                                                      &
@@ -703,8 +703,4 @@
              IF ( nest_domain )  THEN
                 CALL pmci_ensure_nest_mass_conservation
-!
-!--             pmc_update_new is not necessary if nesting is made at each
-!--             substep
-                CALL pmci_update_new
              ENDIF
 
