Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1796)
+++ palm/trunk/SOURCE/Makefile	(revision 1797)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# dependcy for check_for_restart updated
 # 
 # Former revisions:
@@ -342,5 +342,5 @@
 calc_precipitation.o: modules.o mod_kinds.o
 calc_radiation.o: modules.o mod_kinds.o
-check_for_restart.o: modules.o mod_kinds.o
+check_for_restart.o: modules.o mod_kinds.o pmc_interface.o
 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o netcdf_interface.o
 check_parameters.o: modules.o mod_kinds.o land_surface_model.o \
Index: palm/trunk/SOURCE/check_for_restart.f90
===================================================================
--- palm/trunk/SOURCE/check_for_restart.f90	(revision 1796)
+++ palm/trunk/SOURCE/check_for_restart.f90	(revision 1797)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! check now accounts for nesting mode
 ! 
 ! Former revisions:
@@ -70,13 +70,19 @@
                termination_time_needed, time_restart,                          &
                time_since_reference_point, write_binary
+
     USE kinds
+
     USE pegrid
 
+    USE pmc_interface,                                                         &
+        ONLY:  comm_world_nesting, cpl_id, nested_run
+
     IMPLICIT NONE
 
-
-    LOGICAL :: terminate_run_l           !<
-    LOGICAL :: do_stop_now = .FALSE.     !<
-    LOGICAL :: do_restart_now = .FALSE.  !<
+    INTEGER ::  global_communicator       !< global communicator to be used here
+
+    LOGICAL ::  terminate_run_l           !<
+    LOGICAL ::  do_stop_now = .FALSE.     !<
+    LOGICAL ::  do_restart_now = .FALSE.  !<
 
     REAL(wp) ::  remaining_time !<
@@ -96,11 +102,20 @@
     ENDIF
 
+!
+!-- Set the global communicator to be used (depends on the mode in which PALM is
+!-- running)
+    IF ( nested_run )  THEN
+       global_communicator = comm_world_nesting
+    ELSE
+       global_communicator = comm2d
+    ENDIF
+
 #if defined( __parallel )
 !
 !-- Make a logical OR for all processes. Stop the model run if at least
-!-- one processor has reached the time limit.
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,        &
-                        MPI_LOR, comm2d, ierr )
+!-- one process has reached the time limit.
+    IF ( collective_wait )  CALL MPI_BARRIER( global_communicator, ierr )
+    CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,     &
+                        MPI_LOR, global_communicator, ierr )
 #else
     terminate_run = terminate_run_l
@@ -164,9 +179,9 @@
 #if defined( __parallel )
 !
-!-- Make a logical OR for all processes. Stop the model run if at least
-!-- one processor has reached the time limit.
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+!-- Make a logical OR for all processes. Stop the model run if a flag file has
+!-- been detected above.
+    IF ( collective_wait )  CALL MPI_BARRIER( global_communicator, ierr )
     CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL,        &
-                        MPI_LOR, comm2d, ierr )
+                        MPI_LOR, global_communicator, ierr )
 #else
     terminate_run = terminate_run_l
@@ -256,5 +271,6 @@
 !-- the start of a continuation run, except if the user forced to stop the
 !-- run without restart
-    IF ( terminate_run  .AND.  myid == 0 .AND. .NOT. do_stop_now)  THEN
+    IF ( terminate_run  .AND.  myid == 0  .AND.  cpl_id == 1  .AND.            &
+         .NOT. do_stop_now)  THEN
 
        OPEN ( 90, FILE='CONTINUE_RUN', FORM='FORMATTED' )
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1796)
+++ palm/trunk/SOURCE/header.f90	(revision 1797)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-!
+! output of nesting datatransfer mode
 !
 ! Former revisions:
@@ -291,5 +291,5 @@
 
     USE pmc_interface,                                                         &
-        ONLY:  nested_run, nesting_mode
+        ONLY:  nested_run, nesting_datatransfer_mode, nesting_mode
 
     USE radiation_model_mod,                                                   &
@@ -483,5 +483,6 @@
     IF ( nested_run )  THEN
 
-       WRITE ( io, 600 )  TRIM( nesting_mode )
+       WRITE ( io, 600 )  TRIM( nesting_mode ),                                &
+                          TRIM( nesting_datatransfer_mode )
        CALL pmc_get_model_info( ncpl = ncpl, cpl_id = my_cpl_id )
 
@@ -2476,5 +2477,6 @@
 600 FORMAT (/' Nesting informations:'/ &
             ' --------------------'/ &
-            ' Nesting mode:                     ',A// &
+            ' Nesting mode:                     ',A/ &
+            ' Nesting-datatransfer mode:        ',A// &
             ' Nest id  parent  number   lower left coordinates   name'/ &
             ' (*=me)     id    of PEs      x (m)     y (m)' )
Index: palm/trunk/SOURCE/pmc_client.f90
===================================================================
--- palm/trunk/SOURCE/pmc_client.f90	(revision 1796)
+++ palm/trunk/SOURCE/pmc_client.f90	(revision 1797)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! introduction of different datatransfer modes
 !
 ! Former revisions:
@@ -401,5 +401,5 @@
        do i=1,me%inter_npes
           aPE => me%PEs(i)
-          ar  => aPE%array_list(next_array_in_list)    !actual array is last array in list
+          ar  => aPE%array_list(next_array_in_list)
           ar%NrDims    = NrDims
           ar%A_dim     = dims
@@ -575,8 +575,15 @@
       INTEGER(kind=MPI_ADDRESS_KIND)    ::  target_disp
 
-      t1 = PMC_Time()
-      CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for server to fill buffer
-      t2 = PMC_Time()-t1
-      if(present(WaitTime)) WaitTime = t2
+!
+!--   Synchronization of the model is done in pmci_client_synchronize and pmci_server_synchronize
+!--   Therefor the RMA window can be filled without sychronization at this point and a barrier
+!--   is not necessary
+!--   Please note that WaitTime has to be set in PMC_S_FillBuffer AND PMC_C_GetBuffer
+      if(present(WaitTime))  then
+         t1 = PMC_Time()
+         CALL MPI_Barrier(me%intra_comm, ierr)
+         t2 = PMC_Time()
+         WaitTime = t2-t1
+      end if
 
       CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for buffer is filled
Index: palm/trunk/SOURCE/pmc_handle_communicator.f90
===================================================================
--- palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1796)
+++ palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1797)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! introduction of different datatransfer modes,
+! export of comm_world_nesting
 !
 ! Former revisions:
@@ -90,4 +91,5 @@
    ! Coupler Setup
 
+   INTEGER                                    :: m_world_comm !global nesting communicator
    INTEGER                                    :: m_my_CPL_id  !Coupler id of this model
    INTEGER                                    :: m_Parent_id  !Coupler id of parent of this model
@@ -125,5 +127,6 @@
  CONTAINS
 
-   SUBROUTINE pmc_init_model( comm, nesting_mode, pmc_status )
+   SUBROUTINE pmc_init_model( comm, nesting_datatransfer_mode, nesting_mode,   &
+                              pmc_status )
 
       USE control_parameters,                                                  &
@@ -136,4 +139,5 @@
 
       CHARACTER(LEN=7), INTENT(OUT) ::  nesting_mode
+      CHARACTER(LEN=7), INTENT(OUT) ::  nesting_datatransfer_mode
 
       INTEGER, INTENT(OUT)                ::  comm
@@ -148,4 +152,5 @@
       pmc_status   = pmc_status_ok
       comm         = -1
+      m_world_comm = MPI_COMM_WORLD
       m_my_cpl_id  = -1
       clientcount  =  0
@@ -159,5 +164,6 @@
       IF ( m_world_rank == 0 )  THEN
 
-         CALL read_coupling_layout( nesting_mode, pmc_status )
+         CALL read_coupling_layout( nesting_datatransfer_mode, nesting_mode,   &
+                                    pmc_status )
 
          IF ( pmc_status /= pmc_no_namelist_found  .AND.                       &
@@ -221,5 +227,7 @@
          CALL MPI_BCAST( m_couplers(i)%lower_left_x, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
          CALL MPI_BCAST( m_couplers(i)%lower_left_y, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
+!--      TO_DO: the next two calls can and should probably moved outside this loop
          CALL MPI_BCAST( nesting_mode, LEN( nesting_mode ), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
+         CALL MPI_BCAST( nesting_datatransfer_mode, LEN(nesting_datatransfer_mode), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
       ENDDO
 
@@ -329,7 +337,7 @@
 !
 !-- 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 )
+    SUBROUTINE pmc_get_model_info( comm_world_nesting, cpl_id, cpl_name,       &
+                                   cpl_parent_id, lower_left_x, lower_left_y,  &
+                                   ncpl, npe_total, request_for_cpl_id )
 
       USE kinds
@@ -341,4 +349,5 @@
       INTEGER, INTENT(IN), OPTIONAL ::  request_for_cpl_id
 
+      INTEGER, INTENT(OUT), OPTIONAL ::  comm_world_nesting
       INTEGER, INTENT(OUT), OPTIONAL ::  cpl_id
       INTEGER, INTENT(OUT), OPTIONAL ::  cpl_parent_id
@@ -364,6 +373,9 @@
 !
 !--   Return the requested information
+      IF ( PRESENT( comm_world_nesting )  )  THEN
+         comm_world_nesting = m_world_comm
+      ENDIF
       IF ( PRESENT( cpl_id )        )  THEN
-         cpl_id        = requested_cpl_id
+         cpl_id = requested_cpl_id
       ENDIF
       IF ( PRESENT( cpl_parent_id ) )  THEN
@@ -371,17 +383,17 @@
       ENDIF
       IF ( PRESENT( cpl_name )      )  THEN
-         cpl_name      = m_couplers(requested_cpl_id)%name
+         cpl_name = m_couplers(requested_cpl_id)%name
       ENDIF
       IF ( PRESENT( ncpl )          )  THEN
-         ncpl          = m_ncpl
+         ncpl = m_ncpl
       ENDIF
       IF ( PRESENT( npe_total )     )  THEN
-         npe_total     = m_couplers(requested_cpl_id)%npe_total
+         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
+         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
+         lower_left_y = m_couplers(requested_cpl_id)%lower_left_y
       ENDIF
 
@@ -400,9 +412,11 @@
 
 
- SUBROUTINE read_coupling_layout( nesting_mode, pmc_status )
+ SUBROUTINE read_coupling_layout( nesting_datatransfer_mode, nesting_mode,     &
+                                  pmc_status )
 
     IMPLICIT NONE
 
-    CHARACTER(LEN=7) ::  nesting_mode
+    CHARACTER(LEN=7), INTENT(INOUT) ::  nesting_mode
+    CHARACTER(LEN=7), INTENT(INOUT) ::  nesting_datatransfer_mode
 
     INTEGER, INTENT(INOUT) ::  pmc_status
@@ -412,5 +426,5 @@
 
 
-    NAMELIST /nestpar/  domain_layouts, nesting_mode
+    NAMELIST /nestpar/  domain_layouts, nesting_datatransfer_mode, nesting_mode
 
 !
Index: palm/trunk/SOURCE/pmc_interface.f90
===================================================================
--- palm/trunk/SOURCE/pmc_interface.f90	(revision 1796)
+++ palm/trunk/SOURCE/pmc_interface.f90	(revision 1797)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! ------------------
-!
+! introduction of different datatransfer modes,
+! further formatting cleanup, parameter checks added (including mismatches
+! between root and client model settings),
+! +routine pmci_check_setting_mismatches
+! comm_world_nesting introduced
 !
 ! Former revisions:
@@ -145,4 +149,5 @@
 !
 !-- Coupler setup
+    INTEGER(iwp), SAVE      ::  comm_world_nesting     !:
     INTEGER(iwp), SAVE      ::  cpl_id  = 1            !:
     CHARACTER(LEN=32), SAVE ::  cpl_name               !:
@@ -152,4 +157,7 @@
 !
 !-- Control parameters, will be made input parameters later
+    CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !: steering
+                                                         !: parameter for data-
+                                                         !: transfer mode
     CHARACTER(LEN=7), SAVE ::  nesting_mode = 'two-way'  !: steering parameter
                                                          !: for 1- or 2-way nesting
@@ -308,6 +316,6 @@
 
 
-    INTERFACE pmci_client_datatrans
-       MODULE PROCEDURE pmci_client_datatrans
+    INTERFACE pmci_check_setting_mismatches
+       MODULE PROCEDURE pmci_check_setting_mismatches
     END INTERFACE
 
@@ -320,4 +328,8 @@
     END INTERFACE
 
+    INTERFACE pmci_datatrans
+       MODULE PROCEDURE pmci_datatrans
+    END INTERFACE pmci_datatrans
+
     INTERFACE pmci_ensure_nest_mass_conservation
        MODULE PROCEDURE pmci_ensure_nest_mass_conservation
@@ -346,13 +358,13 @@
     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
+           anterp_relax_length_t, client_to_server, comm_world_nesting,        &
+           cpl_id, nested_run, nesting_datatransfer_mode, nesting_mode,        &
+           server_to_client
     PUBLIC pmci_client_initialize
     PUBLIC pmci_client_synchronize
+    PUBLIC pmci_datatrans
     PUBLIC pmci_ensure_nest_mass_conservation
     PUBLIC pmci_init
     PUBLIC pmci_modelconfiguration
-    PUBLIC pmci_server_datatrans
     PUBLIC pmci_server_initialize
     PUBLIC pmci_server_synchronize
@@ -365,4 +377,7 @@
  SUBROUTINE pmci_init( world_comm )
 
+    USE control_parameters,                                                  &
+        ONLY:  message_string
+
     IMPLICIT NONE
 
@@ -376,5 +391,6 @@
 
 
-    CALL pmc_init_model( world_comm, nesting_mode, pmc_status )
+    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,  &
+                         pmc_status )
 
     IF ( pmc_status == pmc_no_namelist_found )  THEN
@@ -390,4 +406,22 @@
 
 !
+!-- Check steering parameter values
+    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                              &
+         TRIM( nesting_mode ) /= 'two-way' )                                   &
+    THEN
+       message_string = 'illegal nesting mode: ' // TRIM( nesting_mode )
+       CALL message( 'pmci_init', 'PA0417', 3, 2, 0, 6, 0 )
+    ENDIF
+
+    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                 &
+         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                 &
+         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                      &
+    THEN
+       message_string = 'illegal nesting datatransfer mode: '                  &
+                        // TRIM( nesting_datatransfer_mode )
+       CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
+    ENDIF
+
+!
 !-- Set the general steering switch which tells PALM that its a nested run
     nested_run = .TRUE.
@@ -396,5 +430,6 @@
 !-- 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_model_info( cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,   &
+    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,          &
+                             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,                &
@@ -445,4 +480,8 @@
 !-- Initialize PMC Server
     CALL pmci_setup_server
+!
+!-- Check for mismatches between settings of master and client variables
+!-- (e.g., all clients have to follow the end_time settings of the root master)
+    CALL pmci_check_setting_mismatches
 
     CALL location_message( 'finished', .TRUE. )
@@ -701,5 +740,5 @@
        ELSE
 !
-!--       TO_DO: Klaus: comment why thie dummy allocation is required
+!--       TO_DO: Klaus: comment why this dummy allocation is required
           ALLOCATE( index_list(6,1) )
           CALL pmc_s_set_2d_index_list( client_id, index_list )
@@ -789,8 +828,4 @@
                                      SIZE(define_coarse_grid_real), 0, 21, ierr )
           CALL pmc_recv_from_server( define_coarse_grid_int,  3, 0, 22, ierr )
-
-!
-!--       Receive also the dz-,zu- and zw-arrays here.
-!--       TO_DO: what is the meaning of above comment
 !
 !--        Debug-printouts - keep them
@@ -852,6 +887,8 @@
        
 !
-!--    TO_DO: give comments what is happening here
+!--    Find the index bounds for the nest domain in the coarse-grid index space
        CALL pmci_map_fine_to_coarse_grid
+!
+!--    TO_DO: Klaus give a comment what is happening here
        CALL pmc_c_get_2d_index_list
 
@@ -879,6 +916,7 @@
 
 !
-!--    Two-way coupling
-!--    TO_DO: comment what is happening here
+!--    Two-way coupling.
+!--    Precompute the index arrays and relaxation functions for the
+!--    anterpolation
        IF ( nesting_mode == 'two-way' )  THEN
           CALL pmci_init_anterp_tophat
@@ -895,5 +933,7 @@
 
     SUBROUTINE pmci_map_fine_to_coarse_grid
-
+!
+!--    Determine index bounds of interpolation/anterpolation area in the coarse
+!--    grid index space
        IMPLICIT NONE
 
@@ -907,6 +947,5 @@
 
 !
-!--    Determine indices of interpolation/anterpolation area in the coarse grid
-!--    If the fine- and coarse grid nodes do not match.
+!--    If the fine- and coarse grid nodes do not match:
        loffset = MOD( coord_x(nxl), cg%dx )
        xexl    = cg%dx + loffset
@@ -991,5 +1030,5 @@
 !--    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.
+!--    and interp_tril_t.
 
        IMPLICIT NONE
@@ -1199,4 +1238,5 @@
 !--    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) )
@@ -1235,4 +1275,5 @@
 
           ENDDO
+
        ENDIF
 
@@ -1240,4 +1281,5 @@
 !--    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) )
@@ -1252,5 +1294,5 @@
           DO  j = nys, nyn
 !
-!--          Right boundary for u.
+!--          Right boundary for u
              i   = nxr + 1
              kb  = nzb_u_inner(j,i)
@@ -1262,7 +1304,6 @@
              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.
+!
+!--          Right boundary for v
              i   = nxr + 1
              kb  = nzb_v_inner(j,i)
@@ -1276,4 +1317,5 @@
 
           ENDDO
+
        ENDIF
 
@@ -1281,4 +1323,5 @@
 !--    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) )
@@ -1291,7 +1334,8 @@
           direction      = 1
           inc            = 1
+
           DO  i = nxl, nxr
 !
-!--          South boundary for u.
+!--          South boundary for u
              j   = -1
              kb  =  nzb_u_inner(j,i)
@@ -1303,5 +1347,4 @@
              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
@@ -1315,5 +1358,7 @@
              logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
              lcr(0:ncorr-1) = 1.0_wp
-          ENDDO
+
+          ENDDO
+
        ENDIF
 
@@ -1321,4 +1366,5 @@
 !--    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) )
@@ -1331,7 +1377,8 @@
           direction      = 1
           inc            = 1
+
           DO  i = nxl, nxr
 !
-!--          North boundary for u.
+!--          North boundary for u
              j   = nyn + 1
              kb  = nzb_u_inner(j,i)
@@ -1343,7 +1390,6 @@
              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.
+!
+!--          North boundary for v
              j   = nyn + 1
              kb  = nzb_v_inner(j,i)
@@ -1357,13 +1403,16 @@
 
           ENDDO
+
        ENDIF
 
 !       
-!--    Then vertical walls and corners if necessary.
+!--    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,       &
@@ -1374,5 +1423,4 @@
              DO  j = nys, nyn
                 DO  k = nzb, nzt_topo_nestbc_l
-
 !
 !--                Wall for u on the south side, but not on the north side
@@ -1392,285 +1440,317 @@
                       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
-                     inc        = -1                 
-                     wall_index =  j + 1
-                     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 w on the south side, but not on the north side.
-                  i  = -1
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.    &
-                       ( nzb_w_outer(j,i) == nzb_w_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_w_l(k,j,2) = inc * lc
-                     logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the north side, but not on the south side.
-                  i  = -1
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND.     &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
-                     inc        = -1
-                     wall_index =  j+1 
-                     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_w_l(k,j,2) = inc * lc
-                     logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF   !  IF ( nest_bound_l )
+                   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 + 1
+                      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 w on the south side, but not on the north side.
+                   i  = -1
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_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_w_l(k,j,2) = inc * lc
+                      logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the north side, but not on the south side.
+                   i  = -1
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
+                      inc        = -1
+                      wall_index =  j+1
+                      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_w_l(k,j,2) = inc * lc
+                      logc_ratio_w_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+          ENDIF   !  IF ( nest_bound_l )
 
 !        
-!--      Right boundary.
-         IF ( nest_bound_r )  THEN
-            ALLOCATE( logc_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
-            ALLOCATE( logc_ratio_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
-            logc_w_r       = 0
-            logc_ratio_w_r = 1.0_wp
-            direction      = 2
-            i  = nxr + 1        
-            DO  j = nys, nyn
-               DO  k = nzb, nzt_topo_nestbc_r
-
-!
-!--               Wall for u on the south side, but not on the north side.
-                  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_r(k,j,2) = inc * lc
-                     logc_ratio_u_r(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.
-                  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+1
-                     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_r(k,j,2) = inc * lc
-                     logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the south side, but not on the north side.
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.    &
-                       ( nzb_w_outer(j,i) == nzb_w_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_w_r(k,j,2) = inc * lc
-                     logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-!
-!--               Wall for w on the north side, but not on the south side.
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) .AND.     &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
-                     inc        = -1
-                     wall_index =  j+1 
-                     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_w_r(k,j,2) = inc * lc
-                     logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF   !  IF ( nest_bound_r )
+!--       Right boundary
+          IF ( nest_bound_r )  THEN
+
+             ALLOCATE( logc_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
+             ALLOCATE( logc_ratio_w_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,       &
+                                      0:ncorr-1) )
+             logc_w_r       = 0
+             logc_ratio_w_r = 1.0_wp
+             direction      = 2
+             i  = nxr + 1
+
+             DO  j = nys, nyn
+                DO  k = nzb, nzt_topo_nestbc_r
+!
+!--                Wall for u on the south side, but not on the north side
+                   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_r(k,j,2) = inc * lc
+                      logc_ratio_u_r(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
+                   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+1
+                      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_r(k,j,2) = inc * lc
+                      logc_ratio_u_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the south side, but not on the north side
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_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_w_r(k,j,2) = inc * lc
+                      logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the north side, but not on the south side
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j+1,i) ) )  THEN
+                      inc        = -1
+                      wall_index =  j+1
+                      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_w_r(k,j,2) = inc * lc
+                      logc_ratio_w_r(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+          ENDIF   !  IF ( nest_bound_r )
 
 !        
-!--      South boundary.
-         IF ( nest_bound_s )  THEN
-            ALLOCATE( logc_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2) )
-            ALLOCATE( logc_ratio_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2, 0:ncorr-1) )
-            logc_w_s       = 0
-            logc_ratio_w_s = 1.0_wp  
-            direction      = 3
-            DO  i = nxl, nxr
-               DO  k = nzb, nzt_topo_nestbc_s
-
-!
-!--               Wall for v on the left side, but not on the right side.
-                  j  = 0
-                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.    &
-                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
-                     inc        =  1
-                     wall_index =  i 
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_v_s(k,i,2) = inc * lc
-                     logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-!
-!--               Wall for v on the right side, but not on the left side.
-                  j  = 0
-                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.    &
-                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
-                     inc        = -1
-                     wall_index =  i+1
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_v_s(k,i,2) = inc * lc
-                     logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the left side, but not on the right side.
-                  j  = -1
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.    &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
-                     inc        =  1
-                     wall_index =  i
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_w_s(k,i,2) = inc * lc
-                     logc_ratio_w_s(k,i,2,0:ncorr - 1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the right side, but not on the left side.
-                  j  = -1
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.    &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
-                     inc        = -1
-                     wall_index =  i+1 
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_w_s(k,i,2) = inc * lc
-                     logc_ratio_w_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF   !  IF (nest_bound_s )
+!--       South boundary
+          IF ( nest_bound_s )  THEN
+
+             ALLOCATE( logc_w_s(nzb:nzt_topo_nestbc_s, nxl:nxr, 1:2) )
+             ALLOCATE( logc_ratio_w_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,       &
+                                      0:ncorr-1) )
+             logc_w_s       = 0
+             logc_ratio_w_s = 1.0_wp
+             direction      = 3
+
+             DO  i = nxl, nxr
+                DO  k = nzb, nzt_topo_nestbc_s
+!
+!--                Wall for v on the left side, but not on the right side
+                   j  = 0
+                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.       &
+                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
+                      inc        =  1
+                      wall_index =  i
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_v_s(k,i,2) = inc * lc
+                      logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for v on the right side, but not on the left side
+                   j  = 0
+                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.       &
+                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
+                      inc        = -1
+                      wall_index =  i+1
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_v_s(k,i,2) = inc * lc
+                      logc_ratio_v_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the left side, but not on the right side
+                   j  = -1
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
+                      inc        =  1
+                      wall_index =  i
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_w_s(k,i,2) = inc * lc
+                      logc_ratio_w_s(k,i,2,0:ncorr - 1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the right side, but not on the left side
+                   j  = -1
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
+                      inc        = -1
+                      wall_index =  i+1
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_w_s(k,i,2) = inc * lc
+                      logc_ratio_w_s(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+          ENDIF   !  IF (nest_bound_s )
 
 !        
-!--      North boundary.
-         IF ( nest_bound_n )  THEN
-            ALLOCATE( logc_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2) )
-            ALLOCATE( logc_ratio_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2, 0:ncorr-1) )
-            logc_w_n       = 0
-            logc_ratio_w_n = 1.0_wp
-            direction      = 3
-            j  = nyn + 1
-            DO  i = nxl, nxr
-               DO  k = nzb, nzt_topo_nestbc_n
-
-!
-!--               Wall for v on the left side, but not on the right side.
-                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.    &
-                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
-                     inc        = 1
-                     wall_index = i 
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_v_n(k,i,2) = inc * lc
-                     logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for v on the right side, but not on the left side.
-                  IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.    &
-                       ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
-                     inc        = -1                 
-                     wall_index =  i + 1
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_v_n(k,i,2) = inc * lc
-                     logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the left side, but not on the right side.
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.    &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
-                     inc        = 1
-                     wall_index = i
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_w_n(k,i,2) = inc * lc
-                     logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-
-!
-!--               Wall for w on the right side, but not on the left side.
-                  IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) .AND.     &
-                       ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
-                     inc        = -1
-                     wall_index =  i+1 
-                     CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, i, 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_w_n(k,i,2) = inc * lc
-                     logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
-                     lcr(0:ncorr-1) = 1.0_wp
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF   !  IF ( nest_bound_n )
-      ENDIF   !  IF ( topography /= 'flat' )
-
-   END SUBROUTINE pmci_init_loglaw_correction
+!--       North boundary
+          IF ( nest_bound_n )  THEN
+
+             ALLOCATE( logc_w_n(nzb:nzt_topo_nestbc_n, nxl:nxr, 1:2) )
+             ALLOCATE( logc_ratio_w_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,       &
+                                      0:ncorr-1) )
+             logc_w_n       = 0
+             logc_ratio_w_n = 1.0_wp
+             direction      = 3
+             j  = nyn + 1
+
+             DO  i = nxl, nxr
+                DO  k = nzb, nzt_topo_nestbc_n
+!
+!--                Wall for v on the left side, but not on the right side
+                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) )  .AND.       &
+                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i-1) ) )  THEN
+                      inc        = 1
+                      wall_index = i
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_v_n(k,i,2) = inc * lc
+                      logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for v on the right side, but not on the left side
+                   IF ( ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) )  .AND.       &
+                        ( nzb_v_outer(j,i) == nzb_v_outer(j,i+1) ) )  THEN
+                      inc        = -1
+                      wall_index =  i + 1
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_v_n(k,i,2) = inc * lc
+                      logc_ratio_v_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the left side, but not on the right side
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i-1) ) )  THEN
+                      inc        = 1
+                      wall_index = i
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_w_n(k,i,2) = inc * lc
+                      logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+!
+!--                Wall for w on the right side, but not on the left side
+                   IF ( ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) )  .AND.       &
+                        ( nzb_w_outer(j,i) == nzb_w_outer(j,i+1) ) )  THEN
+                      inc        = -1
+                      wall_index =  i+1
+                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
+                          k, i, 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_w_n(k,i,2) = inc * lc
+                      logc_ratio_w_n(k,i,2,0:ncorr-1) = lcr(0:ncorr-1)
+                      lcr(0:ncorr-1) = 1.0_wp
+                   ENDIF
+
+                ENDDO
+             ENDDO
+
+          ENDIF   !  IF ( nest_bound_n )
+
+       ENDIF   !  IF ( topography /= 'flat' )
+
+    END SUBROUTINE pmci_init_loglaw_correction
 
 
@@ -1901,321 +1981,337 @@
 
 
-!-- TO_DO:  indentation and wrap long lines from here on to the end of the file
-   SUBROUTINE pmci_init_anterp_tophat
-!
-!--   Precomputation of the client-array indices for 
-!--   corresponding coarse-grid array index and the 
-!--   Under-relaxation coefficients to be used by anterp_tophat.
-
-      IMPLICIT NONE
-
-      INTEGER(iwp) ::  i        !: Fine-grid index
-      INTEGER(iwp) ::  ii       !: Coarse-grid index
-      INTEGER(iwp) ::  istart   !:
-      INTEGER(iwp) ::  j        !: Fine-grid index
-      INTEGER(iwp) ::  jj       !: Coarse-grid index
-      INTEGER(iwp) ::  jstart   !:
-      INTEGER(iwp) ::  k        !: Fine-grid index
-      INTEGER(iwp) ::  kk       !: Coarse-grid index
-      INTEGER(iwp) ::  kstart   !:
-      REAL(wp)     ::  xi       !:
-      REAL(wp)     ::  eta      !:
-      REAL(wp)     ::  zeta     !:
+
+
+    SUBROUTINE pmci_init_anterp_tophat
+!
+!--    Precomputation of the client-array indices for
+!--    corresponding coarse-grid array index and the
+!--    Under-relaxation coefficients to be used by anterp_tophat.
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i        !: Fine-grid index
+       INTEGER(iwp) ::  ii       !: Coarse-grid index
+       INTEGER(iwp) ::  istart   !:
+       INTEGER(iwp) ::  j        !: Fine-grid index
+       INTEGER(iwp) ::  jj       !: Coarse-grid index
+       INTEGER(iwp) ::  jstart   !:
+       INTEGER(iwp) ::  k        !: Fine-grid index
+       INTEGER(iwp) ::  kk       !: Coarse-grid index
+       INTEGER(iwp) ::  kstart   !:
+       REAL(wp)     ::  xi       !:
+       REAL(wp)     ::  eta      !:
+       REAL(wp)     ::  zeta     !:
      
-
-!
-!--   Default values:
-      IF ( anterp_relax_length_l < 0.0_wp )  THEN
-         anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx
-      ENDIF
-      IF ( anterp_relax_length_r < 0.0_wp )  THEN
-         anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx
-      ENDIF
-      IF ( anterp_relax_length_s < 0.0_wp )  THEN
-         anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy
-      ENDIF
-      IF ( anterp_relax_length_n < 0.0_wp )  THEN
-         anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy
-      ENDIF
-      IF ( anterp_relax_length_t < 0.0_wp )  THEN
-         anterp_relax_length_t = 0.1_wp * zu(nzt)
-      ENDIF
-
-!
-!--   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
-      kctu = kk - 1 
-
-      kk = 0
-      DO  WHILE ( cg%zw(kk) < zw(nzt-1) )
-         kk = kk + 1
-      ENDDO
-      kctw = kk - 1 
-
-      ALLOCATE( iflu(icl:icr) )
-      ALLOCATE( iflo(icl:icr) )
-      ALLOCATE( ifuu(icl:icr) )
-      ALLOCATE( ifuo(icl:icr) )
-      ALLOCATE( jflv(jcs:jcn) )
-      ALLOCATE( jflo(jcs:jcn) )
-      ALLOCATE( jfuv(jcs:jcn) )
-      ALLOCATE( jfuo(jcs:jcn) )
-      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  ii = icl, icr  
-         i = istart
-         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
-            i = i + 1
-         ENDDO
-         iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
-         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) ) 
-            i = i + 1
-         ENDDO
-         ifuu(ii) = MIN( MAX( i, nxlg ), nxrg )
-         istart = iflu(ii)
-      ENDDO
-
-!
-!--   i-indices of others for each l-index value.
-      istart = nxlg
-      DO  ii = icl, icr   
-         i = istart
-         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(ii) )  .AND.  ( i < nxrg ) ) 
-            i = i + 1
-         ENDDO
-         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(ii) = MIN(MAX(i,nxlg),nxrg)
-         istart = iflo(ii)
-      ENDDO
-
-!
-!--   j-indices of v for each m-index value.
-      jstart = nysg
-      DO  jj = jcs, jcn
-         j = jstart
-         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
-            j = j + 1
-         ENDDO
-         jflv(jj) = MIN( MAX( j, nysg ), nyng ) 
-         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) ) 
-            j = j + 1
-         ENDDO
-         jfuv(jj) = MIN( MAX( j, nysg ), nyng )
-         jstart = jflv(jj)
-      ENDDO
-
-!
-!--   j-indices of others for each m-index value.
-      jstart = nysg
-      DO  jj = jcs, jcn
-         j = jstart
-         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(jj) )  .AND.  ( j < nyng ) )  
-            j = j + 1
-         ENDDO
-         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(jj) = MIN( MAX( j, nysg ), nyng )
-         jstart = jflv(jj)
-      ENDDO
-
-!
-!--   k-indices of w for each n-index value.
-      kstart  = 0
-      kflw(0) = 0
-      kfuw(0) = 0
-      DO  kk = 1, kctw
-         k = kstart
-         DO  WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) )  .AND.  ( k < nzt ) ) 
-            k = k + 1
-         ENDDO
-         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(kk) = MIN( MAX( k, 1 ), nzt + 1 )
-         kstart = kflw(kk)
-      ENDDO
-
-!
-!--   k-indices of others for each n-index value.
-      kstart  = 0
-      kflo(0) = 0
-      kfuo(0) = 0
-      DO  kk = 1, kctu
-         k = kstart
-         DO  WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) )  .AND.  ( k < nzt ) ) 
-            k = k + 1
-         ENDDO
-         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(kk) = MIN( MAX( k-1, 1 ), nzt + 1 )
-         kstart = kflo(kk)
-      ENDDO
+!
+!--    Default values:
+       IF ( anterp_relax_length_l < 0.0_wp )  THEN
+          anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx
+       ENDIF
+       IF ( anterp_relax_length_r < 0.0_wp )  THEN
+          anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx
+       ENDIF
+       IF ( anterp_relax_length_s < 0.0_wp )  THEN
+          anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy
+       ENDIF
+       IF ( anterp_relax_length_n < 0.0_wp )  THEN
+          anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy
+       ENDIF
+       IF ( anterp_relax_length_t < 0.0_wp )  THEN
+          anterp_relax_length_t = 0.1_wp * zu(nzt)
+       ENDIF
+
+!
+!--    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
+       kctu = kk - 1
+
+       kk = 0
+       DO  WHILE ( cg%zw(kk) < zw(nzt-1) )
+          kk = kk + 1
+       ENDDO
+       kctw = kk - 1
+
+       ALLOCATE( iflu(icl:icr) )
+       ALLOCATE( iflo(icl:icr) )
+       ALLOCATE( ifuu(icl:icr) )
+       ALLOCATE( ifuo(icl:icr) )
+       ALLOCATE( jflv(jcs:jcn) )
+       ALLOCATE( jflo(jcs:jcn) )
+       ALLOCATE( jfuv(jcs:jcn) )
+       ALLOCATE( jfuo(jcs:jcn) )
+       ALLOCATE( kflw(0:kctw) )
+       ALLOCATE( kflo(0:kctu) )
+       ALLOCATE( kfuw(0:kctw) )
+       ALLOCATE( kfuo(0:kctu) )
+
+!
+!--    i-indices of u for each ii-index value
+       istart = nxlg
+       DO  ii = icl, icr
+          i = istart
+          DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  &
+                      ( i < nxrg ) )
+             i = i + 1
+          ENDDO
+          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
+          DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.  &
+                      ( i < nxrg ) )
+             i = i + 1
+          ENDDO
+          ifuu(ii) = MIN( MAX( i, nxlg ), nxrg )
+          istart = iflu(ii)
+       ENDDO
+
+!
+!--    i-indices of others for each ii-index value
+       istart = nxlg
+       DO  ii = icl, icr
+          i = istart
+          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) )  .AND.     &
+                      ( i < nxrg ) )
+             i = i + 1
+          ENDDO
+          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(ii) = MIN(MAX(i,nxlg),nxrg)
+          istart = iflo(ii)
+       ENDDO
+
+!
+!--    j-indices of v for each jj-index value
+       jstart = nysg
+       DO  jj = jcs, jcn
+          j = jstart
+          DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  &
+                      ( j < nyng ) )
+             j = j + 1
+          ENDDO
+          jflv(jj) = MIN( MAX( j, nysg ), nyng )
+          DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.  &
+                      ( j < nyng ) )
+             j = j + 1
+          ENDDO
+          jfuv(jj) = MIN( MAX( j, nysg ), nyng )
+          jstart = jflv(jj)
+       ENDDO
+
+!
+!--    j-indices of others for each jj-index value
+       jstart = nysg
+       DO  jj = jcs, jcn
+          j = jstart
+          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) )  .AND.     &
+                      ( j < nyng ) )
+             j = j + 1
+          ENDDO
+          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(jj) = MIN( MAX( j, nysg ), nyng )
+          jstart = jflv(jj)
+       ENDDO
+
+!
+!--    k-indices of w for each kk-index value
+       kstart  = 0
+       kflw(0) = 0
+       kfuw(0) = 0
+       DO  kk = 1, kctw
+          k = kstart
+          DO  WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) )  .AND.       &
+                      ( k < nzt ) )
+             k = k + 1
+          ENDDO
+          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(kk) = MIN( MAX( k, 1 ), nzt + 1 )
+          kstart = kflw(kk)
+       ENDDO
+
+!
+!--    k-indices of others for each kk-index value
+       kstart  = 0
+       kflo(0) = 0
+       kfuo(0) = 0
+       DO  kk = 1, kctu
+          k = kstart
+          DO  WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) )  .AND.       &
+                      ( k < nzt ) )
+             k = k + 1
+          ENDDO
+          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(kk) = MIN( MAX( k-1, 1 ), nzt + 1 )
+          kstart = kflo(kk)
+       ENDDO
      
 !
-!--   Spatial under-relaxation coefficients
-      ALLOCATE( frax(icl:icr) )
-
-      DO  ii = icl, icr
-         IF ( nest_bound_l ) THEN
-            xi   = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) / anterp_relax_length_l )**4
-         ELSEIF ( nest_bound_r ) THEN
-            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(ii)  = xi / ( 1.0_wp + xi )
-      ENDDO
-
-      ALLOCATE( fray(jcs:jcn) )
-
-      DO  jj = jcs, jcn
-         IF ( nest_bound_s ) THEN           
-            eta  = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) / anterp_relax_length_s )**4
-         ELSEIF ( nest_bound_n ) THEN
-            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(jj)  = eta / ( 1.0_wp + eta )
-      ENDDO
+!--    Spatial under-relaxation coefficients
+       ALLOCATE( frax(icl:icr) )
+
+       DO  ii = icl, icr
+          IF ( nest_bound_l )  THEN
+             xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) /   &
+                    anterp_relax_length_l )**4
+          ELSEIF ( nest_bound_r )  THEN
+             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(ii) = xi / ( 1.0_wp + xi )
+       ENDDO
+
+       ALLOCATE( fray(jcs:jcn) )
+
+       DO  jj = jcs, jcn
+          IF ( nest_bound_s )  THEN
+             eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) /  &
+                     anterp_relax_length_s )**4
+          ELSEIF ( nest_bound_n )  THEN
+             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(jj) = eta / ( 1.0_wp + eta )
+       ENDDO
       
-      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
-
-   END SUBROUTINE pmci_init_anterp_tophat
-
-
-
-   SUBROUTINE pmci_init_tkefactor
-
-!
-!--   Computes the scaling factor for the SGS TKE from coarse grid to be used 
-!--   as BC for the fine grid. Based on the Kolmogorov energy spectrum
-!--   for the inertial subrange and assumption of sharp cut-off of the resolved 
-!--   energy spectrum. Near the surface, the reduction of TKE is made 
-!--   smaller than further away from the surface.
-!
-!                      Antti Hellsten 4.3.2015
-!
-!--   Extended for non-flat topography and variable dz.
-! 
-!                      Antti Hellsten 26.3.2015
-!
-!--   The current near-wall adaptation can be replaced by a new one which
-!--   uses a step function [0,1] based on the logc-arrays. AH 30.12.2015
-      IMPLICIT NONE
-      REAL(wp), PARAMETER  ::  cfw = 0.2_wp          !:
-      REAL(wp), PARAMETER  ::  c_tkef = 0.6_wp       !:
-      REAL(wp)             ::  fw                    !:
-      REAL(wp), PARAMETER  ::  fw0 = 0.9_wp          !:
-      REAL(wp)             ::  glsf                  !:
-      REAL(wp)             ::  glsc                  !:
-      REAL(wp)             ::  height                !:
-      REAL(wp), PARAMETER  ::  p13 = 1.0_wp/3.0_wp   !:
-      REAL(wp), PARAMETER  ::  p23 = 2.0_wp/3.0_wp   !:
-      INTEGER(iwp)         ::  k                     !:
-      INTEGER(iwp)         ::  kc                    !:
+       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
+
+    END SUBROUTINE pmci_init_anterp_tophat
+
+
+
+    SUBROUTINE pmci_init_tkefactor
+
+!
+!--    Computes the scaling factor for the SGS TKE from coarse grid to be used
+!--    as BC for the fine grid. Based on the Kolmogorov energy spectrum
+!--    for the inertial subrange and assumption of sharp cut-off of the resolved
+!--    energy spectrum. Near the surface, the reduction of TKE is made
+!--    smaller than further away from the surface.
+
+       IMPLICIT NONE
+       REAL(wp), PARAMETER ::  cfw = 0.2_wp          !:
+       REAL(wp), PARAMETER ::  c_tkef = 0.6_wp       !:
+       REAL(wp)            ::  fw                    !:
+       REAL(wp), PARAMETER ::  fw0 = 0.9_wp          !:
+       REAL(wp)            ::  glsf                  !:
+       REAL(wp)            ::  glsc                  !:
+       REAL(wp)            ::  height                !:
+       REAL(wp), PARAMETER ::  p13 = 1.0_wp/3.0_wp   !:
+       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !:
+       INTEGER(iwp)        ::  k                     !:
+       INTEGER(iwp)        ::  kc                    !:
         
 
-      IF ( nest_bound_l )  THEN
-         ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )    
-         tkefactor_l = 0.0_wp
-         i = nxl - 1 
-         DO  j = nysg, nyng
-            DO  k = nzb_s_inner(j,i) + 1, nzt
-               kc     = kco(k+1)
-               glsf   = ( dx * dy * dzu(k) )**p13
-               glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
-               height = zu(k) - zu(nzb_s_inner(j,i))
-               fw     = EXP( -cfw * height / glsf )
-               tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
-            ENDDO
-            tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0
-         ENDDO
-      ENDIF
-
-      IF ( nest_bound_r )  THEN
-         ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )    
-         tkefactor_r = 0.0_wp
-         i = nxr + 1
-         DO  j = nysg, nyng
-            DO  k = nzb_s_inner(j,i) + 1, nzt
-               kc     = kco(k+1)
-               glsf   = ( dx * dy * dzu(k) )**p13
-               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-               height = zu(k) - zu(nzb_s_inner(j,i))
-               fw     = EXP( -cfw * height / glsf )
-               tkefactor_r(k,j) = c_tkef * (fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
-            ENDDO
-            tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0
-         ENDDO
-      ENDIF
+       IF ( nest_bound_l )  THEN
+          ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
+          tkefactor_l = 0.0_wp
+          i = nxl - 1
+          DO  j = nysg, nyng
+             DO  k = nzb_s_inner(j,i) + 1, nzt
+                kc     = kco(k+1)
+                glsf   = ( dx * dy * dzu(k) )**p13
+                glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
+                height = zu(k) - zu(nzb_s_inner(j,i))
+                fw     = EXP( -cfw * height / glsf )
+                tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                                              ( glsf / glsc )**p23 )
+             ENDDO
+             tkefactor_l(nzb_s_inner(j,i),j) = c_tkef * fw0
+          ENDDO
+       ENDIF
+
+       IF ( nest_bound_r )  THEN
+          ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
+          tkefactor_r = 0.0_wp
+          i = nxr + 1
+          DO  j = nysg, nyng
+             DO  k = nzb_s_inner(j,i) + 1, nzt
+                kc     = kco(k+1)
+                glsf   = ( dx * dy * dzu(k) )**p13
+                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
+                height = zu(k) - zu(nzb_s_inner(j,i))
+                fw     = EXP( -cfw * height / glsf )
+                tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                                              ( glsf / glsc )**p23 )
+             ENDDO
+             tkefactor_r(nzb_s_inner(j,i),j) = c_tkef * fw0
+          ENDDO
+       ENDIF
 
       IF ( nest_bound_s )  THEN
-         ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
-         tkefactor_s = 0.0_wp
-         j = nys - 1
-         DO  i = nxlg, nxrg
-            DO  k = nzb_s_inner(j,i) + 1, nzt
-               kc     = kco(k+1)
-               glsf   = ( dx * dy * dzu(k) )**p13
-               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
-               height = zu(k) - zu(nzb_s_inner(j,i))
-               fw     = EXP( -cfw*height / glsf )
-               tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
-            ENDDO
-            tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0
-         ENDDO
-      ENDIF
-
-      IF ( nest_bound_n )  THEN
-         ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
-         tkefactor_n = 0.0_wp    
-         j = nyn + 1
-         DO  i = nxlg, nxrg
-            DO  k = nzb_s_inner(j,i)+1, nzt
-               kc     = kco(k+1)
-               glsf   = ( dx * dy * dzu(k) )**p13
-               glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-               height = zu(k) - zu(nzb_s_inner(j,i))
-               fw     = EXP( -cfw * height / glsf )
-               tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
-            ENDDO
-            tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0
-         ENDDO
-      ENDIF
-
-      ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
-      k = nzt
-      DO  i = nxlg, nxrg
-         DO  j = nysg, nyng
-            kc     = kco(k+1)
-            glsf   = ( dx * dy * dzu(k) )**p13
-            glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
-            height = zu(k) - zu(nzb_s_inner(j,i))
-            fw     = EXP( -cfw * height / glsf )
-            tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) * ( glsf / glsc )**p23 )
-         ENDDO
-      ENDDO
+          ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
+          tkefactor_s = 0.0_wp
+          j = nys - 1
+          DO  i = nxlg, nxrg
+             DO  k = nzb_s_inner(j,i) + 1, nzt
+                kc     = kco(k+1)
+                glsf   = ( dx * dy * dzu(k) )**p13
+                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
+                height = zu(k) - zu(nzb_s_inner(j,i))
+                fw     = EXP( -cfw*height / glsf )
+                tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                                              ( glsf / glsc )**p23 )
+             ENDDO
+             tkefactor_s(nzb_s_inner(j,i),i) = c_tkef * fw0
+          ENDDO
+       ENDIF
+
+       IF ( nest_bound_n )  THEN
+          ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
+          tkefactor_n = 0.0_wp
+          j = nyn + 1
+          DO  i = nxlg, nxrg
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                kc     = kco(k+1)
+                glsf   = ( dx * dy * dzu(k) )**p13
+                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
+                height = zu(k) - zu(nzb_s_inner(j,i))
+                fw     = EXP( -cfw * height / glsf )
+                tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
+                                              ( glsf / glsc )**p23 )
+             ENDDO
+             tkefactor_n(nzb_s_inner(j,i),i) = c_tkef * fw0
+          ENDDO
+       ENDIF
+
+       ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
+       k = nzt
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             kc     = kco(k+1)
+             glsf   = ( dx * dy * dzu(k) )**p13
+             glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
+             height = zu(k) - zu(nzb_s_inner(j,i))
+             fw     = EXP( -cfw * height / glsf )
+             tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
+                                           ( glsf / glsc )**p23 )
+          ENDDO
+       ENDDO
       
     END SUBROUTINE pmci_init_tkefactor
@@ -2366,7 +2462,5 @@
 
 !
-!-- 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.
+!-- List of array names, which can be coupled
     IF ( TRIM( name ) == "u" )  THEN
        IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
@@ -2406,5 +2500,5 @@
        ELSE
 !
-!--       Avoid others to continue
+!--       Prevent others from continuing
           CALL MPI_BARRIER( comm2d, ierr )
        ENDIF
@@ -2650,4 +2744,102 @@
 
 
+ SUBROUTINE pmci_check_setting_mismatches
+!
+!-- Check for mismatches between settings of master and client variables
+!-- (e.g., all clients have to follow the end_time settings of the root model).
+!-- The root model overwrites variables in the other models, so these variables
+!-- only need to be set once in file PARIN.
+
+#if defined( __parallel )
+
+    USE control_parameters,                                                    &
+        ONLY:  dt_restart, end_time, message_string, restart_time, time_restart
+
+    IMPLICIT NONE
+
+    INTEGER ::  ierr
+
+    REAL(wp) ::  dt_restart_root
+    REAL(wp) ::  end_time_root
+    REAL(wp) ::  restart_time_root
+    REAL(wp) ::  time_restart_root
+
+!
+!-- Check the time to be simulated.
+!-- Here, and in the following, the root process communicates the respective
+!-- variable to all others, and its value will then be compared with the local
+!-- values.
+    IF ( pmc_is_rootmodel() )  end_time_root = end_time
+    CALL MPI_BCAST( end_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
+
+    IF ( .NOT. pmc_is_rootmodel() )  THEN
+       IF ( end_time /= end_time_root )  THEN
+          WRITE( message_string, * )  'mismatch between root model and ',      &
+               'client settings &   end_time(root) = ', end_time_root,         &
+               ' &   end_time(client) = ', end_time, ' & client value is set', &
+               ' to root value'
+          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
+                        0 )
+          end_time = end_time_root
+       ENDIF
+    ENDIF
+
+!
+!-- Same for restart time
+    IF ( pmc_is_rootmodel() )  restart_time_root = restart_time
+    CALL MPI_BCAST( restart_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
+
+    IF ( .NOT. pmc_is_rootmodel() )  THEN
+       IF ( restart_time /= restart_time_root )  THEN
+          WRITE( message_string, * )  'mismatch between root model and ',      &
+               'client settings &   restart_time(root) = ', restart_time_root, &
+               ' &   restart_time(client) = ', restart_time, ' & client ',     &
+               'value is set to root value'
+          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
+                        0 )
+          restart_time = restart_time_root
+       ENDIF
+    ENDIF
+
+!
+!-- Same for dt_restart
+    IF ( pmc_is_rootmodel() )  dt_restart_root = dt_restart
+    CALL MPI_BCAST( dt_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
+
+    IF ( .NOT. pmc_is_rootmodel() )  THEN
+       IF ( dt_restart /= dt_restart_root )  THEN
+          WRITE( message_string, * )  'mismatch between root model and ',      &
+               'client settings &   dt_restart(root) = ', dt_restart_root,     &
+               ' &   dt_restart(client) = ', dt_restart, ' & client ',         &
+               'value is set to root value'
+          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
+                        0 )
+          dt_restart = dt_restart_root
+       ENDIF
+    ENDIF
+
+!
+!-- Same for time_restart
+    IF ( pmc_is_rootmodel() )  time_restart_root = time_restart
+    CALL MPI_BCAST( time_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
+
+    IF ( .NOT. pmc_is_rootmodel() )  THEN
+       IF ( time_restart /= time_restart_root )  THEN
+          WRITE( message_string, * )  'mismatch between root model and ',      &
+               'client settings &   time_restart(root) = ', time_restart_root, &
+               ' &   time_restart(client) = ', time_restart, ' & client ',     &
+               'value is set to root value'
+          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
+                        0 )
+          time_restart = time_restart_root
+       ENDIF
+    ENDIF
+
+#endif
+
+ END SUBROUTINE pmci_check_setting_mismatches
+
+
+
  SUBROUTINE pmci_ensure_nest_mass_conservation
 
@@ -2787,4 +2979,7 @@
     REAL(wp), DIMENSION(1) ::  dtl         !:
 
+
+    CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )
+
 !
 !-- First find the smallest native time step of all the clients of the current
@@ -2813,4 +3008,6 @@
     ENDDO
 
+    CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )
+
 #endif
  END SUBROUTINE pmci_server_synchronize
@@ -2836,4 +3033,7 @@
     dtl(1) = dt_3d
     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 )
@@ -2845,4 +3045,7 @@
 !--    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 sync', 'stop' )
+
     ENDIF
 
@@ -2853,4 +3056,7 @@
 
  SUBROUTINE pmci_set_swaplevel( swaplevel )
+!
+!-- After each Runge-Kutta sub-timestep, alternately set buffer one or buffer
+!-- two active
 
     IMPLICIT NONE
@@ -2862,6 +3068,4 @@
     INTEGER(iwp)            ::  m          !:
 
-!
-!-- After each timestep, alternately set buffer one or buffer two active
     DO  m = 1, SIZE( pmc_server_for_client )-1
        client_id = pmc_server_for_client(m)
@@ -2870,4 +3074,56 @@
 
  END SUBROUTINE pmci_set_swaplevel
+
+
+
+ SUBROUTINE pmci_datatrans( local_nesting_mode )
+!
+!-- Althoug nesting_mode is a variable of this model, pass it as an argument to
+!-- allow for example to force one-way initialization phase
+
+    IMPLICIT NONE
+
+    INTEGER(iwp)           ::  ierr   !:
+    INTEGER(iwp)           ::  istat  !:
+
+    CHARACTER(LEN=*),INTENT(IN) ::  local_nesting_mode
+
+    IF ( local_nesting_mode == 'one-way' )  THEN
+
+       CALL pmci_client_datatrans( server_to_client )
+       CALL pmci_server_datatrans( server_to_client )
+
+    ELSE
+
+       IF( nesting_datatransfer_mode == 'cascade' )  THEN
+
+          CALL pmci_client_datatrans( server_to_client )
+          CALL pmci_server_datatrans( server_to_client )
+
+          CALL pmci_server_datatrans( client_to_server )
+          CALL pmci_client_datatrans( client_to_server )
+
+       ELSEIF( nesting_datatransfer_mode == 'overlap')  THEN
+
+          CALL pmci_server_datatrans( server_to_client )
+          CALL pmci_client_datatrans( server_to_client )
+
+          CALL pmci_client_datatrans( client_to_server )
+          CALL pmci_server_datatrans( client_to_server )
+
+       ELSEIF( TRIM( nesting_datatransfer_mode ) == 'mixed' )  THEN
+
+          CALL pmci_server_datatrans( server_to_client )
+          CALL pmci_client_datatrans( server_to_client )
+
+          CALL pmci_server_datatrans( client_to_server )
+          CALL pmci_client_datatrans( client_to_server )
+
+       ENDIF
+
+    ENDIF
+
+ END SUBROUTINE pmci_datatrans
+
 
 
@@ -2890,35 +3146,11 @@
     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)
-       IF ( myid==0 )  THEN
-          CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )
-          dtl(1) = MIN( dtl(1), dtc(1) )
-          dt_3d  = dtl(1)
-       ENDIF
-    ENDDO
-
-!
-!-- 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
        client_id = PMC_Server_for_Client(m)
-       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 sync', 'stop' )
        
        IF ( direction == server_to_client )  THEN
           CALL cpu_log( log_point_s(71), 'pmc server send', 'start' )
-          CALL pmc_s_fillbuffer( client_id, waittime=waittime )
+          CALL pmc_s_fillbuffer( client_id )
           CALL cpu_log( log_point_s(71), 'pmc server send', 'stop' )
        ELSE
@@ -2936,13 +3168,16 @@
 !
 !--          Inside buildings/topography reset velocities and TKE back to zero.
-!--          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
-                   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(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
+                   pt(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
+                   IF ( humidity  .OR.  passive_scalar )  THEN
+                      q(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
+                   ENDIF
                 ENDDO
              ENDDO
@@ -2969,6 +3204,4 @@
     INTEGER(iwp) ::  jcn         !:
     
-    REAL(wp) ::  waittime    !:
-
     REAL(wp), DIMENSION(1) ::  dtl         !:
     REAL(wp), DIMENSION(1) ::  dts         !:
@@ -2977,16 +3210,4 @@
     dtl = dt_3d
     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 )
-          CALL pmc_recv_from_server( dts, SIZE( dts ), 0, 102, ierr )
-          dt_3d = dts(1)
-       ENDIF
-
-!
-!--    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 sync', 'stop' )
-
 !
 !--    Client domain boundaries in the server indice space.
@@ -2999,16 +3220,20 @@
 
           CALL cpu_log( log_point_s(73), 'pmc client recv', 'start' )
-          CALL pmc_c_getbuffer( WaitTime = WaitTime )
+          CALL pmc_c_getbuffer( )
           CALL cpu_log( log_point_s(73), 'pmc client recv', 'stop' )
 
+          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'start' )
           CALL pmci_interpolation
+          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' )
 
        ELSE
 !
-!--       direction == server_to_client
+!--       direction == client_to_server
+          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'start' )
           CALL pmci_anterpolation
+          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'stop' )
 
           CALL cpu_log( log_point_s(74), 'pmc client send', 'start' )
-          CALL pmc_c_putbuffer( WaitTime = WaitTime )
+          CALL pmc_c_putbuffer( )
           CALL cpu_log( log_point_s(74), 'pmc client send', 'stop' )
 
@@ -3234,14 +3459,14 @@
 
       CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
-                               kfuo, nzb_u_inner, 'u' )
+                               kfuo, 'u' )
       CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
-                               kfuo, nzb_v_inner, 'v' )
+                               kfuo, 'v' )
       CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
-                               kfuw, nzb_w_inner, 'w' )
+                               kfuw, 'w' )
       CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo,    &
-                               kfuo, nzb_s_inner, 's' )
+                               kfuo, 's' )
       IF ( humidity  .OR.  passive_scalar )  THEN
          CALL pmci_anterp_tophat( q, qc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
-                                  kfuo, nzb_s_inner, 's' )
+                                  kfuo, 's' )
       ENDIF
 
@@ -3256,26 +3481,28 @@
 !--   Interpolation of ghost-node values used as the client-domain boundary
 !--   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
+!--   based on trilinear interpolation.
+
       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           !:
-      REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2,0:ncorr-1), INTENT(IN) ::  logc_ratio   !: 
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r1x          !:
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r2x          !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r1y          !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r2y          !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r1z          !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r2z          !:
+      REAL(wp), DIMENSION(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(nzb:nzt_topo_nestbc,nys:nyn,1:2,0:ncorr-1),          &
+                                      INTENT(IN)    ::  logc_ratio   !:
+      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r1x     !:
+      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r2x     !:
+      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r1y     !:
+      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r2y     !:
+      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z     !:
+      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z     !:
       
-      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                             ::  ic           !:
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                             ::  jc           !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)                   ::  kb           !:
-      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                             ::  kc           !: 
-      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2), INTENT(IN)       ::  logc         !: 
-      INTEGER(iwp)                                                               ::  nzt_topo_nestbc   !:
+      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic     !:
+      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc     !:
+      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb     !:
+      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !:
+      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nys:nyn,1:2),                &
+                                          INTENT(IN)           ::  logc   !:
+      INTEGER(iwp) ::  nzt_topo_nestbc   !:
 
       CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
@@ -3312,8 +3539,9 @@
       
 ! 
-!--   Check which edge is to be handled: left or right. Note the assumption that the same PE never 
-!--   holds both left and right nest boundaries. Should this be changed?
+!--   Check which edge is to be handled
       IF ( edge == 'l' )  THEN
-         IF ( var == 'u' )  THEN   !  For u, nxl is a ghost node, but not for the other variables.
+!
+!--      For u, nxl is a ghost node, but not for the other variables
+         IF ( var == 'u' )  THEN
             i  = nxl
             ib = nxl - 1 
@@ -3327,6 +3555,6 @@
       ENDIF
       
-      DO  j = nys, nyn + 1
-         DO  k = kb(j,i), nzt + 1
+      DO  j = nys, nyn+1
+         DO  k = kb(j,i), nzt+1
             l = ic(i)
             m = jc(j)
@@ -3345,10 +3573,11 @@
 !--   Generalized log-law-correction algorithm.
 !--   Doubly two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays 
-!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction.
+!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine
+!--   pmci_init_loglaw_correction.
 !
 !--   Solid surface below the node 
       IF ( var == 'u' .OR. var == 'v' )  THEN           
          DO  j = nys, nyn
-            k = kb(j,i) + 1
+            k = kb(j,i)+1
             IF ( ( logc(k,j,1) /= 0 )  .AND.  ( logc(k,j,2) == 0 ) )  THEN
                k1 = logc(k,j,1)
@@ -3362,7 +3591,8 @@
 
 !
-!--   In case of non-flat topography, also vertical walls and corners need to be treated.
-!--   Only single and double wall nodes are corrected. Triple and higher-multiple wall nodes 
-!--   are not corrected as the log law would not be valid anyway in such locations.
+!--   In case of non-flat topography, also vertical walls and corners need to be
+!--   treated. Only single and double wall nodes are corrected. Triple and
+!--   higher-multiple wall nodes are not corrected as the log law would not be
+!--   valid anyway in such locations.
       IF ( topography /= 'flat' )  THEN
          IF ( var == 'u' .OR. var == 'w' )  THEN                 
@@ -3371,12 +3601,13 @@
 !--         Solid surface only on south/north side of the node                   
             DO  j = nys, nyn
-               DO  k = kb(j,i) + 1, nzt_topo_nestbc
+               DO  k = kb(j,i)+1, nzt_topo_nestbc
                   IF ( ( logc(k,j,2) /= 0 )  .AND.  ( logc(k,j,1) == 0 ) )  THEN
 
 !
-!--                  Direction of the wall-normal index is carried in as the sign of logc.            
+!--                  Direction of the wall-normal index is carried in as the
+!--                  sign of logc
                      jinc = SIGN( 1, logc(k,j,2) )
                      j1   = ABS( logc(k,j,2) )
-                     DO  jcorr=0, ncorr - 1
+                     DO  jcorr = 0, ncorr-1
                         jco = j + jinc * jcorr
                         f(k,jco,i) = logc_ratio(k,j,2,jcorr) * f(k,j1,i)
@@ -3396,10 +3627,12 @@
                   jinc = SIGN( 1, logc(k,j,2) )
                   j1   = ABS( logc(k,j,2) )                 
-                  DO  jcorr = 0, ncorr - 1
+                  DO  jcorr = 0, ncorr-1
                      jco = j + jinc * jcorr
-                     DO  kcorr = 0, ncorr - 1
+                     DO  kcorr = 0, ncorr-1
                         kco = k + kcorr
-                        f(kco,jco,i) = 0.5_wp * ( logc_ratio(k,j,1,kcorr) * f(k1,j,i)  &
-                                     + logc_ratio(k,j,2,jcorr) * f(k,j1,i) )
+                        f(kco,jco,i) = 0.5_wp * ( logc_ratio(k,j,1,kcorr) *    &
+                                                  f(k1,j,i)                    &
+                                                + logc_ratio(k,j,2,jcorr) *    &
+                                                  f(k,j1,i) )
                      ENDDO
                   ENDDO
@@ -3420,6 +3653,6 @@
             ENDDO
          ELSEIF ( edge == 'r' )  THEN           
-            DO  j = nys, nyn + 1
-               DO  k = kb(j,i), nzt + 1
+            DO  j = nys, nyn+1
+               DO  k = kb(j,i), nzt+1
                   f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
                ENDDO
@@ -3431,9 +3664,9 @@
 !--   Store the boundary values also into the other redundant ghost node layers
       IF ( edge == 'l' )  THEN
-         DO ibgp = -nbgp, ib
+         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
+      ELSEIF ( edge == 'r' )  THEN
+         DO  ibgp = ib, nx+nbgp
             f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
          ENDDO
@@ -3444,6 +3677,7 @@
 
 
-   SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, &
-                                  nzt_topo_nestbc, edge, var )
+   SUBROUTINE pmci_interp_tril_sn( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
+                                   r2z, kb, logc, logc_ratio,                  &
+                                   nzt_topo_nestbc, edge, var )
 
 !
@@ -3451,43 +3685,27 @@
 !--   conditions. This subroutine handles the south and north 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 horifontal surfaces, not yet near vertical 
-!--   walls.
-!--                                      Antti Hellsten 26.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                 !:
-      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)                 ::  fc                !:
-      REAL(wp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2,0:ncorr-1), INTENT(IN) ::  logc_ratio        !:
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r1x               !:
-      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                                 ::  r2x               !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r1y               !:
-      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                                 ::  r2y               !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r1z               !:
-      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                                 ::  r2z               !:
+
+      REAL(wp), DIMENSION(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(nzb:nzt_topo_nestbc,nxl:nxr,1:2,0:ncorr-1),          &
+                                      INTENT(IN)    ::  logc_ratio    !:
+      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r1x           !:
+      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r2x           !:
+      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r1y           !:
+      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r2y           !:
+      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z           !:
+      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z           !:
       
-      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                             ::  ic                !:
-      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                             ::  jc                !:
-      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)                   ::  kb                !:
-      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                             ::  kc                !:
-      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2), INTENT(IN)       ::  logc              !: 
-      INTEGER(iwp)                                                               ::  nzt_topo_nestbc   !:
+      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !:
+      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !:
+      INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !:
+      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
+      INTEGER(iwp), DIMENSION(nzb:nzt_topo_nestbc,nxl:nxr,1:2),                &
+                                          INTENT(IN)           ::  logc  !:
+      INTEGER(iwp) ::  nzt_topo_nestbc   !:
 
       CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
@@ -3521,8 +3739,9 @@
       
 !
-!--   Check which edge is to be handled: south or north. Note the assumption that the same PE never 
-!--   holds both south and north nest boundaries. Should this be changed?
+!--   Check which edge is to be handled: south or north
       IF ( edge == 's' )  THEN
-         IF ( var == 'v' )  THEN   !  For v, nys is a ghost node, but not for the other variables.
+!
+!--      For v, nys is a ghost node, but not for the other variables
+         IF ( var == 'v' )  THEN
             j  = nys
             jb = nys - 1 
@@ -3536,6 +3755,6 @@
       ENDIF
 
-      DO  i = nxl, nxr + 1
-         DO  k = kb(j,i), nzt + 1
+      DO  i = nxl, nxr+1
+         DO  k = kb(j,i), nzt+1
             l = ic(i)
             m = jc(j)
@@ -3554,5 +3773,6 @@
 !--   Generalized log-law-correction algorithm.
 !--   Multiply two-dimensional index arrays logc(:,:,1:2) and log-ratio arrays 
-!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine pmci_init_loglaw_correction.
+!--   logc_ratio(:,:,1:2,0:ncorr-1) have been precomputed in subroutine
+!--   pmci_init_loglaw_correction.
 !
 !--   Solid surface below the node 
@@ -3562,5 +3782,5 @@
             IF ( ( logc(k,i,1) /= 0 )  .AND.  ( logc(k,i,2) == 0 ) )  THEN
                k1 = logc(k,i,1)
-               DO  kcorr = 0, ncorr - 1
+               DO  kcorr = 0, ncorr-1
                   kco = k + kcorr
                   f(kco,j,i) = logc_ratio(k,i,1,kcorr) * f(k1,j,i)
@@ -3571,8 +3791,9 @@
 
 !
-!--   In case of non-flat topography, also vertical walls and corners need to be treated.
-!--   Only single and double wall nodes are corrected. 
-!--   Triple and higher-multiple wall nodes are not corrected as it would be extremely complicated 
-!--   and the log law would not be valid anyway in such locations.
+!--   In case of non-flat topography, also vertical walls and corners need to be
+!--   treated. Only single and double wall nodes are corrected.
+!--   Triple and higher-multiple wall nodes are not corrected as it would be
+!--   extremely complicated and the log law would not be valid anyway in such
+!--   locations.
       IF ( topography /= 'flat' )  THEN
          IF ( var == 'v' .OR. var == 'w' )  THEN
@@ -3585,8 +3806,9 @@
 
 !
-!--                  Direction of the wall-normal index is carried in as the sign of logc.            
+!--                  Direction of the wall-normal index is carried in as the
+!--                  sign of logc
                      iinc = SIGN( 1, logc(k,i,2) )
                      i1  = ABS( logc(k,i,2) )
-                     DO  icorr = 0, ncorr - 1
+                     DO  icorr = 0, ncorr-1
                         ico = i + iinc * icorr
                         f(k,j,ico) = logc_ratio(k,i,2,icorr) * f(k,j,i1)
@@ -3606,10 +3828,12 @@
                   iinc = SIGN( 1, logc(k,i,2) )
                   i1   = ABS( logc(k,i,2) )
-                  DO  icorr = 0, ncorr - 1
+                  DO  icorr = 0, ncorr-1
                      ico = i + iinc * icorr
-                     DO  kcorr = 0, ncorr - 1
+                     DO  kcorr = 0, ncorr-1
                         kco = k + kcorr
-                        f(kco,i,ico) = 0.5_wp * ( logc_ratio(k,i,1,kcorr) * f(k1,j,i)  &
-                                     + logc_ratio(k,i,2,icorr) * f(k,j,i1) )
+                        f(kco,i,ico) = 0.5_wp * ( logc_ratio(k,i,1,kcorr) *    &
+                                                  f(k1,j,i)  &
+                                                + logc_ratio(k,i,2,icorr) *    &
+                                                  f(k,j,i1) )
                      ENDDO
                   ENDDO
@@ -3625,5 +3849,5 @@
          IF ( edge == 's' )  THEN
             DO  i = nxl, nxr + 1
-               DO  k = kb(j,i), nzt + 1
+               DO  k = kb(j,i), nzt+1
                   f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
                ENDDO
@@ -3631,5 +3855,5 @@
          ELSEIF ( edge == 'n' )  THEN
             DO  i = nxl, nxr + 1
-               DO  k = kb(j,i), nzt + 1
+               DO  k = kb(j,i), nzt+1
                   f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
                ENDDO
@@ -3641,9 +3865,9 @@
 !--   Store the boundary values also into the other redundant ghost node layers
       IF ( edge == 's' )  THEN
-         DO jbgp = -nbgp, jb
+         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
+      ELSEIF ( edge == 'n' )  THEN
+         DO  jbgp = jb, ny+nbgp
             f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
          ENDDO
@@ -3654,5 +3878,6 @@
  
 
-   SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, var )
+   SUBROUTINE pmci_interp_tril_t( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z,  &
+                                  r2z, var )
 
 !
@@ -3660,18 +3885,21 @@
 !--   conditions. This subroutine handles the top boundary. 
 !--   This subroutine is based on trilinear interpolation.
-!--   Constant dz is still assumed.
+
       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   !:
+
+      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(nxlg:nxrg), INTENT(IN) ::  ic    !:
+      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN) ::  jc    !:
+      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN) ::  kc    !:
       
       CHARACTER(LEN=1), INTENT(IN) :: var   !:
@@ -3701,6 +3929,6 @@
       ENDIF
      
-      DO  i = nxl - 1, nxr + 1
-         DO  j = nys - 1, nyn + 1
+      DO  i = nxl-1, nxr+1
+         DO  j = nys-1, nyn+1
             l = ic(i)
             m = jc(j)
@@ -3724,7 +3952,7 @@
 !
 !--   Rescale if f is the TKE.
-!--   It is assumed that the bottom surface never reaches the top 
-!---  boundary of a nest domain. 
-      IF ( var == 'e')  THEN
+!--   It is assumed that the bottom surface never reaches the top boundary of a
+!--   nest domain.
+      IF ( var == 'e' )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
@@ -3746,4 +3974,5 @@
 !--    subroutine handles the left and right boundaries. However, this operation
 !--    is only needed in case of one-way coupling.
+
        IMPLICIT NONE
 
@@ -3785,6 +4014,6 @@
        ENDIF
 
-       DO  j = nys, nyn + 1
-          DO  k = kb(j,i), nzt +1
+       DO  j = nys, nyn+1
+          DO  k = kb(j,i), nzt+1
              vdotnor = outnor * u(k,j,ied)
 !
@@ -3806,5 +4035,5 @@
           ENDDO
        ELSEIF ( edge == 'r' )  THEN
-          DO ibgp = ib, nx + nbgp
+          DO ibgp = ib, nx+nbgp
              f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
           ENDDO
@@ -3822,4 +4051,5 @@
 !--    interpolated values by values extrapolated from the domain. This
 !--    subroutine handles the south and north boundaries.
+
        IMPLICIT NONE
 
@@ -3861,6 +4091,6 @@
        ENDIF
 
-       DO  i = nxl, nxr + 1
-          DO  k = kb(j,i), nzt + 1
+       DO  i = nxl, nxr+1
+          DO  k = kb(j,i), nzt+1
              vdotnor = outnor * v(k,jed,i)
 !
@@ -3878,9 +4108,9 @@
 !--    Store the boundary values also into the redundant ghost node layers.
        IF ( edge == 's' )  THEN
-          DO jbgp = -nbgp, jb
+          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
+          DO  jbgp = jb, ny+nbgp
              f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
           ENDDO
@@ -3896,4 +4126,5 @@
 !--    conditions. This subroutine handles the top boundary. It is based on
 !--    trilinear interpolation.
+
        IMPLICIT NONE
 
@@ -3941,5 +4172,5 @@
 
     SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu,   &
-                                   kb, var )
+                                   var )
 !
 !--    Anterpolation of internal-node values to be used as the server-domain
@@ -3947,4 +4178,5 @@
 !--    integration of the fine-grid values contained within the coarse-grid
 !--    cell.
+
        IMPLICIT NONE
 
@@ -3968,12 +4200,10 @@
        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   !:
+       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(0:kct), INTENT(IN)   ::  kfl   !:
+       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)   ::  kfu   !:
 
 
@@ -4049,6 +4279,8 @@
                    ENDIF
                 ENDIF
-
-                fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +                 &
+!
+!--             TO DO: introduce 3-d coarse grid array for precomputed
+!--             1/REAL(nfc) values
+                fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +               &
                                fra * cellsum / REAL( nfc, KIND = wp )
 
Index: palm/trunk/SOURCE/pmc_server.f90
===================================================================
--- palm/trunk/SOURCE/pmc_server.f90	(revision 1796)
+++ palm/trunk/SOURCE/pmc_server.f90	(revision 1797)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! introduction of different datatransfer modes
 !
 ! Former revisions:
@@ -498,8 +498,15 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  data_3d
 
-      t1 = PMC_Time()
-      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)              ! Wait for buffer empty
-      t2 = PMC_Time()
-      if(present(WaitTime)) WaitTime = t2-t1
+!--   Synchronization of the model is done in pmci_client_synchronize and pmci_server_synchronize
+!--   Therefor the RMA window cann be filled without sychronization at this point and the barrier
+!--   is not necessary
+!--   Please note that WaitTime has to be set in PMC_S_FillBuffer AND PMC_C_GetBuffer
+
+      if(present(WaitTime)) then
+        t1 = PMC_Time()
+        CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)
+        t2 = PMC_Time()
+        WaitTime = t2-t1
+      end if
 
       do ip=1,Clients(ClientId)%inter_npes
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1796)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1797)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! introduction of different datatransfer modes
 !
 ! Former revisions:
@@ -258,6 +258,6 @@
     USE pmc_interface,                                                         &
         ONLY:  client_to_server, nested_run, nesting_mode,                     &
-               pmci_ensure_nest_mass_conservation, pmci_client_datatrans,      &
-               pmci_client_synchronize, pmci_server_datatrans,                 &
+               pmci_client_synchronize, pmci_datatrans,                        &
+               pmci_ensure_nest_mass_conservation,                             &
                pmci_server_synchronize, server_to_client
 
@@ -673,20 +673,15 @@
 
           IF ( nested_run )  THEN
-!
-!--          TO_DO: try to give more meaningful comments here
-!--          Domain nesting
-!--          Note that the nesting operations are omitted intentionally on the
-!--          first two RK-substeps.
+
              CALL cpu_log( log_point(60), 'nesting', 'start' )
 !
-!--          From server to client commmunication ( direction=SERVER_TO_CLIENT )
-             CALL pmci_server_datatrans( server_to_client )
-             CALL pmci_client_datatrans( server_to_client )
+!--          Domain nesting. The data transfer subroutines pmci_server_datatrans
+!--          and pmci_client_datatatrans are called inside the wrapper
+!--          subroutine pmci_datatrans according to the control parameters
+!--          nesting_mode and nesting_datatransfer_mode.
+!--          TO_DO: why is nesting_mode given as a parameter here?
+             CALL pmci_datatrans( nesting_mode )
 
              IF ( nesting_mode == 'two-way' )  THEN
-!
-!--             From client to server commmunication ( direction=CLIENT_TO_SERVER )
-                CALL pmci_server_datatrans( client_to_server )
-                CALL pmci_client_datatrans( client_to_server )
 !
 !--             Exchange_horiz is needed for all server-domains after the
