Changeset 1791 for palm/trunk


Ignore:
Timestamp:
Mar 11, 2016 10:41:25 AM (8 years ago)
Author:
raasch
Message:

output of nesting informations of all domains; filling up redundant ghost points; renaming of variables, etc.; formatting cleanup

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1787 r1791  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# dependencies of header changed
    2323#
    2424# Former revisions:
     
    381381global_min_max.o: modules.o mod_kinds.o
    382382header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model.o\
    383    plant_canopy_model.o pmc_interface.o radiation_model.o spectrum.o \
    384    subsidence.o
     383   plant_canopy_model.o pmc_handle_communicator.o pmc_interface.o \
     384   radiation_model.o spectrum.o subsidence.o
    385385impact_of_latent_heat.o: modules.o mod_kinds.o
    386386inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
  • palm/trunk/SOURCE/header.f90

    r1789 r1791  
    1919! Current revisions:
    2020! -----------------
    21 !
    22 ! 
     21! output of nesting informations of all domains
     22!
    2323! Former revisions:
    2424! -----------------
     
    284284               plant_canopy
    285285
     286    USE pmc_handle_communicator,                                               &
     287        ONLY:  pmc_get_model_info
     288
    286289    USE pmc_interface,                                                         &
    287         ONLY:  cpl_id, cpl_parent_id, cpl_name, lower_left_coord_x,            &
    288                lower_left_coord_y, nested_run, nesting_mode
     290        ONLY:  nested_run, nesting_mode
    289291
    290292    USE radiation_model_mod,                                                   &
     
    311313   
    312314    CHARACTER (LEN=26) ::  ver_rev             !<
     315
     316    CHARACTER (LEN=32) ::  cpl_name            !<
    313317   
    314318    CHARACTER (LEN=40) ::  output_format       !<
     
    338342    CHARACTER (LEN=1), DIMENSION(1:3) ::  dir = (/ 'x', 'y', 'z' /)  !<
    339343
    340     INTEGER(iwp) ::  av        !<
    341     INTEGER(iwp) ::  bh        !<
    342     INTEGER(iwp) ::  blx       !<
    343     INTEGER(iwp) ::  bly       !<
    344     INTEGER(iwp) ::  bxl       !<
    345     INTEGER(iwp) ::  bxr       !<
    346     INTEGER(iwp) ::  byn       !<
    347     INTEGER(iwp) ::  bys       !<
    348     INTEGER(iwp) ::  ch        !<
    349     INTEGER(iwp) ::  count     !<
    350     INTEGER(iwp) ::  cwx       !<
    351     INTEGER(iwp) ::  cwy       !<
    352     INTEGER(iwp) ::  cxl       !<
    353     INTEGER(iwp) ::  cxr       !<
    354     INTEGER(iwp) ::  cyn       !<
    355     INTEGER(iwp) ::  cys       !<
    356     INTEGER(iwp) ::  dim       !<
    357     INTEGER(iwp) ::  i         !<
    358     INTEGER(iwp) ::  io        !<
    359     INTEGER(iwp) ::  j         !<
    360     INTEGER(iwp) ::  k         !<
    361     INTEGER(iwp) ::  l         !<
    362     INTEGER(iwp) ::  ll        !<
    363     INTEGER(iwp) ::  mpi_type  !<
     344    INTEGER(iwp) ::  av             !<
     345    INTEGER(iwp) ::  bh             !<
     346    INTEGER(iwp) ::  blx            !<
     347    INTEGER(iwp) ::  bly            !<
     348    INTEGER(iwp) ::  bxl            !<
     349    INTEGER(iwp) ::  bxr            !<
     350    INTEGER(iwp) ::  byn            !<
     351    INTEGER(iwp) ::  bys            !<
     352    INTEGER(iwp) ::  ch             !<
     353    INTEGER(iwp) ::  count          !<
     354    INTEGER(iwp) ::  cpl_parent_id  !<
     355    INTEGER(iwp) ::  cwx            !<
     356    INTEGER(iwp) ::  cwy            !<
     357    INTEGER(iwp) ::  cxl            !<
     358    INTEGER(iwp) ::  cxr            !<
     359    INTEGER(iwp) ::  cyn            !<
     360    INTEGER(iwp) ::  cys            !<
     361    INTEGER(iwp) ::  dim            !<
     362    INTEGER(iwp) ::  i              !<
     363    INTEGER(iwp) ::  io             !<
     364    INTEGER(iwp) ::  j              !<
     365    INTEGER(iwp) ::  k              !<
     366    INTEGER(iwp) ::  l              !<
     367    INTEGER(iwp) ::  ll             !<
     368    INTEGER(iwp) ::  mpi_type       !<
     369    INTEGER(iwp) ::  my_cpl_id      !<
     370    INTEGER(iwp) ::  n              !<
     371    INTEGER(iwp) ::  ncpl           !<
     372    INTEGER(iwp) ::  npe_total      !<
    364373   
    365374    REAL(wp) ::  canopy_height                    !< canopy height (in m)
    366375    REAL(wp) ::  cpuseconds_per_simulated_second  !<
     376    REAL(wp) ::  lower_left_coord_x               !< x-coordinate of nest domain
     377    REAL(wp) ::  lower_left_coord_y               !< y-coordinate of nest domain
    367378
    368379!
     
    468479!-- Nesting informations
    469480    IF ( nested_run )  THEN
    470        WRITE ( io, 600 )  cpl_id, TRIM( cpl_name ), cpl_parent_id,             &
    471                           nesting_mode, lower_left_coord_x, lower_left_coord_y
     481
     482       WRITE ( io, 600 )  TRIM( nesting_mode )
     483       CALL pmc_get_model_info( ncpl = ncpl, cpl_id = my_cpl_id )
     484
     485       DO  n = 1, ncpl
     486          CALL pmc_get_model_info( request_for_cpl_id = n, cpl_name = cpl_name,&
     487                                   cpl_parent_id = cpl_parent_id,              &
     488                                   lower_left_x = lower_left_coord_x,          &
     489                                   lower_left_y = lower_left_coord_y,          &
     490                                   npe_total = npe_total )
     491          IF ( n == my_cpl_id )  THEN
     492             char1 = '*'
     493          ELSE
     494             char1 = ' '
     495          ENDIF
     496          WRITE ( io, 601 )  TRIM( char1 ), n, cpl_parent_id, npe_total,       &
     497                             lower_left_coord_x, lower_left_coord_y,           &
     498                             TRIM( cpl_name )
     499       ENDDO
    472500    ENDIF
    473501    WRITE ( io, 99 )
     
    24432471513 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order ' // &
    24442472            '+ monotonic adjustment')
    2445 600 FORMAT (/' Nesting informations:'/                                        &
    2446             ' Nest id / name:                   ',I2.2,' / ',A,' (parent id: ',I2.2,')'/ &
    2447             ' Nesting mode:                     ',A/ &
    2448             ' Lower left corner coordinates:    ','x = ',F8.2,' m, y = ',F8.2,' m'/)
     2473600 FORMAT (/' Nesting informations:'/ &
     2474            ' --------------------'/ &
     2475            ' Nesting mode:                     ',A// &
     2476            ' Nest id  parent  number   lower left coordinates   name'/ &
     2477            ' (*=me)     id    of PEs      x (m)     y (m)' )
     2478601 FORMAT (2X,A1,1X,I2.2,6X,I2.2,5X,I5,5X,F8.2,2X,F8.2,5X,A)
    24492479
    24502480 END SUBROUTINE header
  • palm/trunk/SOURCE/pmc_client.f90

    r1787 r1791  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Debug write-statement commented out
    2323!
    2424! Former revisions:
     
    150150        end do
    151151
    152         if(me%model_rank == 0) write(0,'(a,5i6)') 'PMC_ClientInit ',me%model_rank,me%model_npes,me%inter_npes,me%intra_rank
     152!        if(me%model_rank == 0) write(0,'(a,5i6)') 'PMC_ClientInit ',me%model_rank,me%model_npes,me%inter_npes,me%intra_rank
    153153
    154154        return
  • palm/trunk/SOURCE/pmc_handle_communicator.f90

    r1787 r1791  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! m_nrofcpl renamed m_ncpl,
     23! pmc_get_local_model_info renamed pmc_get_model_info, some keywords also
     24! renamed and some added,
     25! debug write-statements commented out
    2326!
    2427! Former revisions:
     
    8689   INTEGER                                    :: m_my_CPL_id  !Coupler id of this model
    8790   INTEGER                                    :: m_Parent_id  !Coupler id of parent of this model
    88    INTEGER                                    :: m_NrOfCpl    !Number of Coupler in layout file
    89    TYPE(PMC_layout),DIMENSION(PMC_MAX_MODELL) :: m_couplers   !Information of all coupler
     91   INTEGER                                    :: m_ncpl       !Number of Couplers in layout file
     92
     93   TYPE(PMC_layout),DIMENSION(PMC_MAX_MODELL) :: m_couplers   !Information of all couplers
    9094
    9195   ! MPI settings
     
    110114   END INTERFACE pmc_is_rootmodel
    111115
    112    INTERFACE PMC_get_local_model_info
    113       MODULE PROCEDURE PMC_get_local_model_info
    114    END INTERFACE PMC_get_local_model_info
    115 
    116    PUBLIC pmc_get_local_model_info, pmc_init_model, pmc_is_rootmodel
     116   INTERFACE pmc_get_model_info
     117      MODULE PROCEDURE pmc_get_model_info
     118   END INTERFACE pmc_get_model_info
     119
     120   PUBLIC pmc_get_model_info, pmc_init_model, pmc_is_rootmodel
    117121
    118122 CONTAINS
     
    160164!--         Calculate start PE of every model
    161165            start_pe(1) = 0
    162             DO  i = 2, m_nrofcpl+1
     166            DO  i = 2, m_ncpl+1
    163167               start_pe(i) = start_pe(i-1) + m_couplers(i-1)%npe_total
    164168            ENDDO
     
    167171!--         The number of cores provided with the run must be the same as the
    168172!--         total sum of cores required by all nest domains
    169             IF ( start_pe(m_nrofcpl+1) /= m_world_npes )  THEN
     173            IF ( start_pe(m_ncpl+1) /= m_world_npes )  THEN
    170174               WRITE ( message_string, '(A,I6,A,I6,A)' )                       &
    171175                               'nesting-setup requires more MPI procs (',      &
    172                                start_pe(m_nrofcpl+1), ') than provided (',     &
     176                               start_pe(m_ncpl+1), ') than provided (',        &
    173177                               m_world_npes,')'
    174178               CALL message( 'pmc_init_model', 'PA0229', 3, 2, 0, 6, 0 )
     
    202206      ENDIF
    203207
    204       CALL MPI_BCAST( m_nrofcpl, 1,          MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
    205       CALL MPI_BCAST( start_pe, m_nrofcpl+1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
     208      CALL MPI_BCAST( m_ncpl,          1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
     209      CALL MPI_BCAST( start_pe, m_ncpl+1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat)
    206210
    207211!
    208212!--   Broadcast coupling layout
    209       DO  i = 1, m_nrofcpl
     213      DO  i = 1, m_ncpl
    210214         CALL MPI_BCAST( m_couplers(i)%name, LEN( m_couplers(i)%name ), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
    211215         CALL MPI_BCAST( m_couplers(i)%id,           1, MPI_INTEGER, 0, MPI_COMM_WORLD, istat )
     
    219223!
    220224!--   Assign global MPI processes to individual models by setting the couple id
    221       DO  i = 1, m_nrofcpl
     225      DO  i = 1, m_ncpl
    222226         IF ( m_world_rank >= start_pe(i)  .AND.  m_world_rank < start_pe(i+1) ) &
    223227         THEN
     
    241245!
    242246!--   Broadcast (from PE 0) the parent id and id of every model
    243       DO  i = 1, m_nrofcpl
     247      DO  i = 1, m_ncpl
    244248         CALL MPI_BCAST( m_couplers(i)%parent_id, 1, MPI_INTEGER, 0,           &
    245249                         MPI_COMM_WORLD, istat )
     
    257261!--   different colors.
    258262!--   The grouping was done above with MPI_COMM_SPLIT
    259       DO  i = 2, m_nrofcpl
     263      DO  i = 2, m_ncpl
    260264
    261265         IF ( m_couplers(i)%parent_id == m_my_cpl_id )  THEN
     
    292296
    293297      clientcount = 0
    294       DO  i = 2, m_nrofcpl
     298      DO  i = 2, m_ncpl
    295299         IF ( activeserver(i) == 1 )  THEN
    296300            clientcount = clientcount + 1
     
    321325
    322326!
    323 !-- Make module private variables available to palm
    324    SUBROUTINE pmc_get_local_model_info( my_cpl_id, my_cpl_parent_id, cpl_name, &
    325                                         npe_total, lower_left_x, lower_left_y )
     327!-- Provide module private variables of the pmc for PALM
     328    SUBROUTINE pmc_get_model_info( cpl_id, cpl_name, cpl_parent_id,            &
     329                                   lower_left_x, lower_left_y, ncpl, npe_total,&
     330                                   request_for_cpl_id )
    326331
    327332      USE kinds
     
    330335
    331336      CHARACTER(LEN=*), INTENT(OUT), OPTIONAL ::  cpl_name
    332       INTEGER, INTENT(OUT), OPTIONAL          ::  my_cpl_id
    333       INTEGER, INTENT(OUT), OPTIONAL          ::  my_cpl_parent_id
    334       INTEGER, INTENT(OUT), OPTIONAL          ::  npe_total
    335       REAL(wp), INTENT(OUT), OPTIONAL         ::  lower_left_x
    336       REAL(wp), INTENT(OUT), OPTIONAL         ::  lower_left_y
    337 
    338       IF ( PRESENT( my_cpl_id )           )  my_cpl_id        = m_my_cpl_id
    339       IF ( PRESENT( my_cpl_parent_id )    )  my_cpl_parent_id = m_couplers(my_cpl_id)%parent_id
    340       IF ( PRESENT( cpl_name )            )  cpl_name         = m_couplers(my_cpl_id)%name
    341       IF ( PRESENT( npe_total )           )  npe_total        = m_couplers(my_cpl_id)%npe_total
    342       IF ( PRESENT( lower_left_x )        )  lower_left_x     = m_couplers(my_cpl_id)%lower_left_x
    343       IF ( PRESENT( lower_left_y )        )  lower_left_y     = m_couplers(my_cpl_id)%lower_left_y
    344 
    345    END SUBROUTINE pmc_get_local_model_info
     337
     338      INTEGER, INTENT(IN), OPTIONAL ::  request_for_cpl_id
     339
     340      INTEGER, INTENT(OUT), OPTIONAL ::  cpl_id
     341      INTEGER, INTENT(OUT), OPTIONAL ::  cpl_parent_id
     342      INTEGER, INTENT(OUT), OPTIONAL ::  ncpl
     343      INTEGER, INTENT(OUT), OPTIONAL ::  npe_total
     344
     345      INTEGER ::  requested_cpl_id
     346
     347      REAL(wp), INTENT(OUT), OPTIONAL ::  lower_left_x
     348      REAL(wp), INTENT(OUT), OPTIONAL ::  lower_left_y
     349
     350!
     351!--   Set the requested coupler id
     352      IF ( PRESENT( request_for_cpl_id ) )  THEN
     353         requested_cpl_id = request_for_cpl_id
     354!
     355!--      Check for allowed range of values
     356         IF ( requested_cpl_id < 1 .OR. requested_cpl_id > m_ncpl )  RETURN
     357      ELSE
     358         requested_cpl_id = m_my_cpl_id
     359      ENDIF
     360
     361!
     362!--   Return the requested information
     363      IF ( PRESENT( cpl_id )        )  THEN
     364         cpl_id        = requested_cpl_id
     365      ENDIF
     366      IF ( PRESENT( cpl_parent_id ) )  THEN
     367         cpl_parent_id = m_couplers(requested_cpl_id)%parent_id
     368      ENDIF
     369      IF ( PRESENT( cpl_name )      )  THEN
     370         cpl_name      = m_couplers(requested_cpl_id)%name
     371      ENDIF
     372      IF ( PRESENT( ncpl )          )  THEN
     373         ncpl          = m_ncpl
     374      ENDIF
     375      IF ( PRESENT( npe_total )     )  THEN
     376         npe_total     = m_couplers(requested_cpl_id)%npe_total
     377      ENDIF
     378      IF ( PRESENT( lower_left_x )  )  THEN
     379         lower_left_x  = m_couplers(requested_cpl_id)%lower_left_x
     380      ENDIF
     381      IF ( PRESENT( lower_left_y )  )  THEN
     382         lower_left_y  = m_couplers(requested_cpl_id)%lower_left_y
     383      ENDIF
     384
     385   END SUBROUTINE pmc_get_model_info
    346386
    347387
     
    364404
    365405    INTEGER, INTENT(INOUT) ::  pmc_status
    366     INTEGER                ::  i, istat, iunit
     406    INTEGER                ::  i, istat
    367407
    368408    TYPE(pmc_layout), DIMENSION(pmc_max_modell) ::  domain_layouts
     
    374414!-- Initialize some coupling variables
    375415    domain_layouts(1:pmc_max_modell)%id = -1
    376     m_nrofcpl =   0
    377     iunit     = 345
     416    m_ncpl =   0
    378417
    379418    pmc_status = pmc_status_ok
     
    412451!-- Get the number of nested models given in the nestpar-NAMELIST
    413452    DO  i = 1, pmc_max_modell
    414 
    415        IF ( m_couplers(i)%id /= -1  .AND.  i <= pmc_max_modell )  THEN
    416           WRITE ( 0, '(A,A,1X,3I7,1X,2F10.2)' )  'Set up Model  ',             &
    417                              TRIM( m_couplers(i)%name ), m_couplers(i)%id,     &
    418                              m_couplers(i)%Parent_id, m_couplers(i)%npe_total, &
    419                              m_couplers(i)%lower_left_x,                       &
    420                              m_couplers(i)%lower_left_y
    421        ELSE
    422 !
    423 !--       When id=-1 is found for the first time, the list of domains is
    424 !--       finished (or latest after pmc_max_modell entries
    425           m_nrofcpl = i - 1
    426           EXIT
     453!
     454!--    When id=-1 is found for the first time, the list of domains is finished
     455       IF ( m_couplers(i)%id == -1  .OR.  i == pmc_max_modell )  THEN
     456          IF ( m_couplers(i)%id == -1 )  THEN
     457             m_ncpl = i - 1
     458             EXIT
     459          ELSE
     460             m_ncpl = pmc_max_modell
     461          ENDIF
    427462       ENDIF
    428463
  • palm/trunk/SOURCE/pmc_interface.f90

    r1784 r1791  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! routine pmci_update_new removed,
     23! pmc_get_local_model_info renamed pmc_get_model_info, some keywords also
     24! renamed,
     25! filling up redundant ghost points introduced,
     26! some index bound variables renamed,
     27! further formatting cleanup
    2328!
    2429! Former revisions:
     
    111116
    112117    USE pmc_handle_communicator,                                               &
    113         ONLY:  pmc_get_local_model_info, pmc_init_model, pmc_is_rootmodel,     &
     118        ONLY:  pmc_get_model_info, pmc_init_model, pmc_is_rootmodel,           &
    114119               pmc_no_namelist_found, pmc_server_for_client
    115120
     
    128133    IMPLICIT NONE
    129134
    130 !-- TO_DO: a lot of lines (including comments) in this file exceed the 80 char
    131 !--        limit. Try to reduce as much as possible
    132 
    133 !-- TO_DO: are all of these variables following now really PUBLIC?
    134 !--        Klaus and I guess they are not
    135     PRIVATE    !:  Note that the default publicity is here set to private.
     135    PRIVATE
    136136
    137137!
    138138!-- Constants
    139     INTEGER(iwp), PARAMETER, PUBLIC ::  client_to_server = 2   !:
    140     INTEGER(iwp), PARAMETER, PUBLIC ::  server_to_client = 1   !:
     139    INTEGER(iwp), PARAMETER ::  client_to_server = 2   !:
     140    INTEGER(iwp), PARAMETER ::  server_to_client = 1   !:
    141141
    142142!
    143143!-- Coupler setup
    144     INTEGER(iwp), PUBLIC, SAVE      ::  cpl_id  = 1            !:
    145     CHARACTER(LEN=32), PUBLIC, SAVE ::  cpl_name               !:
    146     INTEGER(iwp), PUBLIC, SAVE      ::  cpl_npe_total          !:
    147     INTEGER(iwp), PUBLIC, SAVE      ::  cpl_parent_id          !:
     144    INTEGER(iwp), SAVE      ::  cpl_id  = 1            !:
     145    CHARACTER(LEN=32), SAVE ::  cpl_name               !:
     146    INTEGER(iwp), SAVE      ::  cpl_npe_total          !:
     147    INTEGER(iwp), SAVE      ::  cpl_parent_id          !:
    148148
    149149!
    150150!-- Control parameters, will be made input parameters later
    151     CHARACTER(LEN=7), PUBLIC, SAVE ::  nesting_mode = 'two-way'  !: steering parameter for one- or two-way nesting
    152 
    153     LOGICAL, PUBLIC, SAVE ::  nested_run = .FALSE.  !: general switch if nested run or not
    154 
    155     REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_l = -1.0_wp   !:
    156     REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_r = -1.0_wp   !:
    157     REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_s = -1.0_wp   !:
    158     REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_n = -1.0_wp   !:
    159     REAL(wp), PUBLIC, SAVE ::  anterp_relax_length_t = -1.0_wp   !:
     151    CHARACTER(LEN=7), SAVE ::  nesting_mode = 'two-way'  !: steering parameter
     152                                                         !: for 1- or 2-way nesting
     153
     154    LOGICAL, SAVE ::  nested_run = .FALSE.  !: general switch
     155
     156    REAL(wp), SAVE ::  anterp_relax_length_l = -1.0_wp   !:
     157    REAL(wp), SAVE ::  anterp_relax_length_r = -1.0_wp   !:
     158    REAL(wp), SAVE ::  anterp_relax_length_s = -1.0_wp   !:
     159    REAL(wp), SAVE ::  anterp_relax_length_n = -1.0_wp   !:
     160    REAL(wp), SAVE ::  anterp_relax_length_t = -1.0_wp   !:
    160161
    161162!
    162163!-- Geometry
    163     REAL(wp), PUBLIC, SAVE                    ::  area_t               !:
     164    REAL(wp), SAVE                            ::  area_t               !:
    164165    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_x              !:
    165166    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE ::  coord_y              !:
    166     REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_x   !:
    167     REAL(wp), PUBLIC, SAVE                    ::  lower_left_coord_y   !:
     167    REAL(wp), SAVE                            ::  lower_left_coord_x   !:
     168    REAL(wp), SAVE                            ::  lower_left_coord_y   !:
    168169
    169170!
    170171!-- Client coarse data arrays
    171     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ec   !:
    172     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  ptc  !:
    173     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  uc   !:
    174     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  vc   !:
    175     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  wc   !:
    176     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET, PUBLIC ::  qc   !:
    177 
    178     INTEGER(iwp), DIMENSION(5)                          ::  coarse_bound   !:
    179     REAL(wp), PUBLIC, SAVE                              ::  xexl           !:
    180     REAL(wp), PUBLIC, SAVE                              ::  xexr           !:
    181     REAL(wp), PUBLIC, SAVE                              ::  yexs           !:
    182     REAL(wp), PUBLIC, SAVE                              ::  yexn           !:
    183     REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !:
    184     REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !:
    185     REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !:
    186     REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !:
    187     REAL(wp), PUBLIC, SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !:
    188 
    189 !
    190 !-- Client interpolation coefficients and client-array indices to be precomputed and stored.
    191     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !:
    192     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !:
    193     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !:
    194     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !:
    195     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !:
    196     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !:
    197     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !:
    198     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !:
    199     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !:
    200     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !:
    201     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !:
    202     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !:
    203     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !:
    204     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !:
    205     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !:
    206     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !:
    207     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
    208     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
    209 
    210 !
    211 !-- Client index arrays and log-ratio arrays for the log-law near-wall corrections.
    212 !-- These are not truly 3-D arrays but multiply 2-D arrays.
    213     INTEGER(iwp), PUBLIC, SAVE :: ncorr  !: ncorr is the 4th dimension of the log_ratio-arrays
    214     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !:
    215     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !:
    216     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !:
    217     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !:
    218     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !:
    219     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !:
    220     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !:
    221     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !:
    222     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !:
    223     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !:
    224     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !:
    225     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !:
    226     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !:
    227     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !:
    228     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !:
    229     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !:
    230     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !:
    231     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !:
    232     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !:
    233     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !:
    234     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !:
    235     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !:
    236     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
    237     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
     172    INTEGER(iwp), DIMENSION(5)                  ::  coarse_bound   !:
     173
     174    REAL(wp), SAVE                              ::  xexl           !:
     175    REAL(wp), SAVE                              ::  xexr           !:
     176    REAL(wp), SAVE                              ::  yexs           !:
     177    REAL(wp), SAVE                              ::  yexn           !:
     178    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !:
     179    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !:
     180    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !:
     181    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !:
     182    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !:
     183
     184    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec   !:
     185    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc  !:
     186    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uc   !:
     187    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vc   !:
     188    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc   !:
     189    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc   !:
     190
     191!
     192!-- Client interpolation coefficients and client-array indices to be precomputed
     193!-- and stored.
     194    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !:
     195    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !:
     196    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !:
     197    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !:
     198    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !:
     199    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !:
     200    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !:
     201    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !:
     202    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !:
     203    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !:
     204    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !:
     205    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !:
     206    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !:
     207    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !:
     208    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !:
     209    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !:
     210    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !:
     211    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !:
     212
     213!
     214!-- Client index arrays and log-ratio arrays for the log-law near-wall
     215!-- corrections. These are not truly 3-D arrays but multiple 2-D arrays.
     216    INTEGER(iwp), SAVE :: ncorr  !: 4th dimension of the log_ratio-arrays
     217    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !:
     218    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !:
     219    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !:
     220    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !:
     221    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !:
     222    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !:
     223    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !:
     224    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !:
     225    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !:
     226    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !:
     227    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !:
     228    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !:
     229    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !:
     230    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !:
     231    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !:
     232    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !:
     233    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !:
     234    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !:
     235    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !:
     236    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !:
     237    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !:
     238    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !:
     239    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !:
     240    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !:
    238241
    239242!
    240243!-- Upper bounds for k in anterpolation.
    241     INTEGER(iwp), PUBLIC, SAVE ::  kceu   !:
    242     INTEGER(iwp), PUBLIC, SAVE ::  kcew   !:
     244    INTEGER(iwp), SAVE ::  kctu   !:
     245    INTEGER(iwp), SAVE ::  kctw   !:
    243246
    244247!
    245248!-- Upper bound for k in log-law correction in interpolation.
    246     INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_l   !:
    247     INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_n   !:
    248     INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_r   !:
    249     INTEGER(iwp), PUBLIC, SAVE ::  nzt_topo_nestbc_s   !:
     249    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_l   !:
     250    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_n   !:
     251    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_r   !:
     252    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_s   !:
    250253
    251254!
    252255!-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
    253     INTEGER(iwp), PUBLIC, SAVE ::  nhll   !:
    254     INTEGER(iwp), PUBLIC, SAVE ::  nhlr   !:
    255     INTEGER(iwp), PUBLIC, SAVE ::  nhls   !:
    256     INTEGER(iwp), PUBLIC, SAVE ::  nhln   !:
     256    INTEGER(iwp), SAVE ::  nhll   !:
     257    INTEGER(iwp), SAVE ::  nhlr   !:
     258    INTEGER(iwp), SAVE ::  nhls   !:
     259    INTEGER(iwp), SAVE ::  nhln   !:
    257260
    258261!
    259262!-- Spatial under-relaxation coefficients for anterpolation.
    260     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !:
    261     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
    262     REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
     263    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !:
     264    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !:
     265    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !:
    263266
    264267!
    265268!-- Client-array indices to be precomputed and stored for anterpolation.
    266     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !:
    267     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !:
    268     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !:
    269     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !:
    270     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !:
    271     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !:
    272     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !:
    273     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !:
    274     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !:
    275     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !:
    276     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
    277     INTEGER(iwp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
    278 
    279 !
    280 !-- Module private variables.
     269    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !:
     270    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !:
     271    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !:
     272    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !:
     273    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !:
     274    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !:
     275    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !:
     276    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !:
     277    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !:
     278    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !:
     279    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !:
     280    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !:
     281
    281282    INTEGER(iwp), DIMENSION(3)          ::  define_coarse_grid_int    !:
    282283    REAL(wp), DIMENSION(7)              ::  define_coarse_grid_real   !:
     
    301302    END TYPE coarsegrid_def
    302303                                         
    303     TYPE(coarsegrid_def), SAVE             ::  cg   !:
     304    TYPE(coarsegrid_def), SAVE ::  cg   !:
    304305
    305306
     
    340341    END INTERFACE pmci_set_swaplevel
    341342
    342     INTERFACE pmci_update_new
    343        MODULE PROCEDURE pmci_update_new
    344     END INTERFACE
    345 
     343    PUBLIC anterp_relax_length_l, anterp_relax_length_r,                       &
     344           anterp_relax_length_s, anterp_relax_length_n,                       &
     345           anterp_relax_length_t, client_to_server, cpl_id, nested_run,        &
     346           nesting_mode, server_to_client
    346347    PUBLIC pmci_client_datatrans
    347348    PUBLIC pmci_client_initialize
     
    354355    PUBLIC pmci_server_synchronize
    355356    PUBLIC pmci_set_swaplevel
    356     PUBLIC pmci_update_new
    357357
    358358
     
    393393!-- Get some variables required by the pmc-interface (and in some cases in the
    394394!-- PALM code out of the pmci) out of the pmc-core
    395     CALL pmc_get_local_model_info( my_cpl_id = cpl_id,                         &
    396                                    my_cpl_parent_id = cpl_parent_id,           &
    397                                    cpl_name = cpl_name,                        &
    398                                    npe_total = cpl_npe_total,                  &
    399                                    lower_left_x = lower_left_coord_x,          &
    400                                    lower_left_y = lower_left_coord_y )
     395    CALL pmc_get_model_info( cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,   &
     396                             cpl_name = cpl_name, npe_total = cpl_npe_total,   &
     397                             lower_left_x = lower_left_coord_x,                &
     398                             lower_left_y = lower_left_coord_y )
    401399!
    402400!-- Set the steering switch which tells the models that they are nested (of
    403 !-- course the root domain (cpl_id = 1 ) is not nested)
     401!-- course the root domain (cpl_id = 1) is not nested)
    404402    IF ( cpl_id >= 2 )  THEN
    405403       nest_domain = .TRUE.
     
    435433
    436434    CALL location_message( 'setup the nested model configuration', .FALSE. )
    437     CALL pmci_setup_coordinates   !:  Compute absolute coordinates valid for all models
    438     CALL pmci_setup_client        !:  Initialize PMC Client (Must be called before pmc_setup_server)
    439     CALL pmci_setup_server        !:  Initialize PMC Server
     435!
     436!-- Compute absolute coordinates for all models
     437    CALL pmci_setup_coordinates
     438!
     439!-- Initialize the client (must be called before pmc_setup_server)
     440    CALL pmci_setup_client
     441!
     442!-- Initialize PMC Server
     443    CALL pmci_setup_server
     444
    440445    CALL location_message( 'finished', .TRUE. )
    441446
     
    449454    IMPLICIT NONE
    450455
    451     INTEGER(iwp)               ::  client_id        !:
    452     INTEGER(iwp)               ::  ierr             !:
    453     INTEGER(iwp)               ::  i                !:
    454     INTEGER(iwp)               ::  j                !:
    455     INTEGER(iwp)               ::  k                !:
    456     INTEGER(iwp)               ::  m                !:
    457     INTEGER(iwp)               ::  nomatch          !:
    458     INTEGER(iwp)               ::  nx_cl            !:
    459     INTEGER(iwp)               ::  ny_cl            !:
    460     INTEGER(iwp)               ::  nz_cl            !:
    461     INTEGER(iwp), DIMENSION(5) ::  val              !:
    462     REAL(wp), DIMENSION(1)     ::  fval             !:
    463     REAL(wp)                   ::  dx_cl            !:
    464     REAL(wp)                   ::  dy_cl            !:
    465     REAL(wp)                   ::  xez              !:
    466     REAL(wp)                   ::  yez              !:
    467     CHARACTER(len=32)          ::  mychannel
    468     CHARACTER(len=32)          ::  myname
     456    CHARACTER(LEN=32) ::  myname
     457
     458    INTEGER(iwp) ::  client_id        !:
     459    INTEGER(iwp) ::  ierr             !:
     460    INTEGER(iwp) ::  i                !:
     461    INTEGER(iwp) ::  j                !:
     462    INTEGER(iwp) ::  k                !:
     463    INTEGER(iwp) ::  m                !:
     464    INTEGER(iwp) ::  nomatch          !:
     465    INTEGER(iwp) ::  nx_cl            !:
     466    INTEGER(iwp) ::  ny_cl            !:
     467    INTEGER(iwp) ::  nz_cl            !:
     468
     469    INTEGER(iwp), DIMENSION(5) ::  val    !:
     470
     471    REAL(wp) ::  dx_cl            !:
     472    REAL(wp) ::  dy_cl            !:
     473    REAL(wp) ::  xez              !:
     474    REAL(wp) ::  yez              !:
     475
     476    REAL(wp), DIMENSION(1) ::  fval             !:
     477
    469478    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !:
    470479    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !:
     
    472481
    473482!
    474 !   Initialize the PMC server
     483!   Initialize the pmc server
    475484    CALL pmc_serverinit
    476485
     
    478487!-- Get coordinates from all clients
    479488    DO  m = 1, SIZE( pmc_server_for_client ) - 1
     489
    480490       client_id = pmc_server_for_client(m)
    481491       IF ( myid == 0 )  THEN       
     
    510520          CALL pmc_recv_from_client( client_id, cl_coord_y, SIZE( cl_coord_y ),&
    511521                                     0, 12, ierr )
    512           WRITE ( 0, * )  'receive from pmc Client ', client_id, nx_cl, ny_cl
     522!          WRITE ( 0, * )  'receive from pmc Client ', client_id, nx_cl, ny_cl
    513523         
    514524          define_coarse_grid_real(1) = lower_left_coord_x
     
    529539          xez = ( nbgp + 1 ) * dx
    530540          yez = ( nbgp + 1 ) * dy
    531           IF ( cl_coord_x(0) < define_coarse_grid_real(1) + xez )          nomatch = 1
    532           IF ( cl_coord_x(nx_cl + 1) > define_coarse_grid_real(5) - xez )  nomatch = 1
    533           IF ( cl_coord_y(0) < define_coarse_grid_real(2) + yez )          nomatch = 1
    534           IF ( cl_coord_y(ny_cl + 1) > define_coarse_grid_real(6) - yez )  nomatch = 1
     541          IF ( ( cl_coord_x(0) < define_coarse_grid_real(1) + xez )       .OR. &
     542               ( cl_coord_x(nx_cl+1) > define_coarse_grid_real(5) - xez ) .OR. &
     543               ( cl_coord_y(0) < define_coarse_grid_real(2) + yez )       .OR. &
     544               ( cl_coord_y(ny_cl+1) > define_coarse_grid_real(6) - yez ) )    &
     545          THEN
     546             nomatch = 1
     547          ENDIF
    535548
    536549          DEALLOCATE( cl_coord_x )
     
    539552!
    540553!--       Send coarse grid information to client
    541           CALL pmc_send_to_client( client_id, Define_coarse_grid_real,         &
    542                                    SIZE(define_coarse_grid_real), 0,           &
    543                                    21, ierr )
    544           CALL pmc_send_to_client( client_id, Define_coarse_grid_int,  3, 0,   &
     554          CALL pmc_send_to_client( client_id, define_coarse_grid_real,         &
     555                                   SIZE( define_coarse_grid_real ), 0, 21,     &
     556                                   ierr )
     557          CALL pmc_send_to_client( client_id, define_coarse_grid_int,  3, 0,   &
    545558                                   22, ierr )
    546559
    547560!
    548561!--       Send local grid to client
    549           CALL pmc_send_to_client( client_id, coord_x, nx+1+2*nbgp, 0, 24, ierr )
    550           CALL pmc_send_to_client( client_id, coord_y, ny+1+2*nbgp, 0, 25, ierr )
     562          CALL pmc_send_to_client( client_id, coord_x, nx+1+2*nbgp, 0, 24,     &
     563                                   ierr )
     564          CALL pmc_send_to_client( client_id, coord_y, ny+1+2*nbgp, 0, 25,     &
     565                                   ierr )
    551566
    552567!
    553568!--       Also send the dzu-, dzw-, zu- and zw-arrays here
    554           CALL pmc_send_to_client( client_id, dzu, nz_cl + 1, 0, 26, ierr )
    555           CALL pmc_send_to_client( client_id, dzw, nz_cl + 1, 0, 27, ierr )
    556           CALL pmc_send_to_client( client_id, zu,  nz_cl + 2, 0, 28, ierr )
    557           CALL pmc_send_to_client( client_id, zw,  nz_cl + 2, 0, 29, ierr )
     569          CALL pmc_send_to_client( client_id, dzu, nz_cl+1, 0, 26, ierr )
     570          CALL pmc_send_to_client( client_id, dzw, nz_cl+1, 0, 27, ierr )
     571          CALL pmc_send_to_client( client_id, zu,  nz_cl+2, 0, 28, ierr )
     572          CALL pmc_send_to_client( client_id, zw,  nz_cl+2, 0, 29, ierr )
    558573
    559574       ENDIF
     
    567582     
    568583       CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
    569        
     584
     585!
     586!--    TO_DO: Klaus: please give a comment what is done here
    570587       CALL pmci_create_index_list
    571588
    572589!
    573590!--    Include couple arrays into server content
     591!--    TO_DO: Klaus: please give a more meaningful comment
    574592       CALL pmc_s_clear_next_array_list
    575593       DO  WHILE ( pmc_s_getnextarray( client_id, myname ) )
     
    587605       IMPLICIT NONE
    588606
     607       INTEGER(iwp) ::  i                  !:
     608       INTEGER(iwp) ::  ic                 !:
     609       INTEGER(iwp) ::  ierr               !:
     610       INTEGER(iwp) ::  j                  !:
     611       INTEGER(iwp) ::  k                  !:
     612       INTEGER(iwp) ::  m                  !:
     613       INTEGER(iwp) ::  n                  !:
     614       INTEGER(iwp) ::  npx                !:
     615       INTEGER(iwp) ::  npy                !:
     616       INTEGER(iwp) ::  nrx                !:
     617       INTEGER(iwp) ::  nry                !:
     618       INTEGER(iwp) ::  px                 !:
     619       INTEGER(iwp) ::  py                 !:
     620       INTEGER(iwp) ::  server_pe          !:
     621
     622       INTEGER(iwp), DIMENSION(2) ::  scoord             !:
     623       INTEGER(iwp), DIMENSION(2) ::  size_of_array      !:
     624
    589625       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  coarse_bound_all   !:
    590        INTEGER(iwp)                               ::  i                  !:
    591        INTEGER(iwp)                               ::  ic                 !:
    592        INTEGER(iwp)                               ::  ierr               !:
    593626       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list         !:
    594        INTEGER(iwp)                               ::  j                  !:
    595        INTEGER(iwp)                               ::  k                  !:
    596        INTEGER(iwp)                               ::  m                  !:
    597        INTEGER(iwp)                               ::  n                  !:
    598        INTEGER(iwp)                               ::  npx                !:
    599        INTEGER(iwp)                               ::  npy                !:
    600        INTEGER(iwp)                               ::  nrx                !:
    601        INTEGER(iwp)                               ::  nry                !:
    602        INTEGER(iwp)                               ::  px                 !:
    603        INTEGER(iwp)                               ::  py                 !:
    604        INTEGER(iwp), DIMENSION(2)                 ::  scoord             !:
    605        INTEGER(iwp)                               ::  server_pe          !:
    606        INTEGER(iwp), DIMENSION(2)                 ::  size_of_array      !:
    607 
    608627
    609628       IF ( myid == 0 )  THEN
     629!--       TO_DO: Klaus: give more specific comment what size_of_array stands for
    610630          CALL pmc_recv_from_client( client_id, size_of_array, 2, 0, 40, ierr )
    611631          ALLOCATE( coarse_bound_all(size_of_array(1),size_of_array(2)) )
     
    628648          CALL MPI_COMM_SIZE( comm1dx, npx, ierr )
    629649          CALL MPI_COMM_SIZE( comm1dy, npy, ierr )
    630 
    631           nrx = nxr - nxl + 1   !  +1 in index because FORTRAN indexing starts with 1, palm with 0
     650!
     651!--       The +1 in index is because PALM starts with nx=0
     652          nrx = nxr - nxl + 1
    632653          nry = nyn - nys + 1
    633654          ic  = 0
    634           DO  k = 1, size_of_array(2)                                   !  loop over all client PEs
    635              DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)       !  area in y required by actual client PE
    636                 DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)    !  area in x required by actual client PE
     655!
     656!--       Loop over all client PEs
     657          DO  k = 1, size_of_array(2)
     658!
     659!--          Area along y required by actual client PE
     660             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
     661!
     662!--             Area along x required by actual client PE
     663                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
     664
    637665                   px = i / nrx
    638666                   py = j / nry
     
    642670                 
    643671                   ic = ic + 1
    644                    index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp       !  First index in Server Array
    645                    index_list(2,ic) = j - ( py * nry ) + 1 + nbgp       !  Second index in Server Array
    646                    index_list(3,ic) = i - coarse_bound_all(1,k) + 1     !  x Index client coarse grid
    647                    index_list(4,ic) = j - coarse_bound_all(3,k) + 1     !  y Index client coarse grid
    648                    index_list(5,ic) = k - 1                             !  PE Number client
    649                    index_list(6,ic) = server_pe                         !  PE Number server
     672!
     673!--                First index in server array
     674                   index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp
     675!
     676!--                Second index in server array
     677                   index_list(2,ic) = j - ( py * nry ) + 1 + nbgp
     678!
     679!--                x index of client coarse grid
     680                   index_list(3,ic) = i - coarse_bound_all(1,k) + 1
     681!
     682!--                y index of client coarse grid
     683                   index_list(4,ic) = j - coarse_bound_all(3,k) + 1
     684!
     685!--                PE number of client
     686                   index_list(5,ic) = k - 1
     687!
     688!--                PE number of server
     689                   index_list(6,ic) = server_pe
     690
    650691                ENDDO
    651692             ENDDO
    652693          ENDDO
     694!
     695!--       TO_DO: Klaus: comment what is done here
    653696          CALL pmc_s_set_2d_index_list( client_id, index_list(:,1:ic) )
     697
    654698       ELSE
    655           ALLOCATE( index_list(6,1) )    !  Dummy allocate
     699!
     700!--       TO_DO: Klaus: comment why thie dummy allocation is required
     701          ALLOCATE( index_list(6,1) )
    656702          CALL pmc_s_set_2d_index_list( client_id, index_list )
    657703       ENDIF
     
    671717    IMPLICIT NONE
    672718
    673     CHARACTER(LEN=DA_Namelen)  ::  myname     !:
    674 
    675     INTEGER(iwp)               ::  i          !:
    676     INTEGER(iwp)               ::  ierr       !:
    677     INTEGER(iwp)               ::  icl        !:
    678     INTEGER(iwp)               ::  icr        !:
    679     INTEGER(iwp)               ::  j          !:
    680     INTEGER(iwp)               ::  jcn        !:
    681     INTEGER(iwp)               ::  jcs        !:
     719    CHARACTER(LEN=da_namelen) ::  myname     !:
     720
     721    INTEGER(iwp) ::  i          !:
     722    INTEGER(iwp) ::  ierr       !:
     723    INTEGER(iwp) ::  icl        !:
     724    INTEGER(iwp) ::  icr        !:
     725    INTEGER(iwp) ::  j          !:
     726    INTEGER(iwp) ::  jcn        !:
     727    INTEGER(iwp) ::  jcs        !:
     728
    682729    INTEGER(iwp), DIMENSION(5) ::  val        !:
    683730   
    684     REAL(wp), DIMENSION(1)     ::  fval       !:
    685     REAL(wp)                   ::  xcs        !:
    686     REAL(wp)                   ::  xce        !:
    687     REAL(wp)                   ::  ycs        !:
    688     REAL(wp)                   ::  yce        !:
     731    REAL(wp) ::  xcs        !:
     732    REAL(wp) ::  xce        !:
     733    REAL(wp) ::  ycs        !:
     734    REAL(wp) ::  yce        !:
     735
     736    REAL(wp), DIMENSION(1) ::  fval       !:
    689737                                             
    690 
     738!
    691739!-- TO_DO: describe what is happening in this if-clause
    692740!-- Root Model does not have Server and is not a client
    693741    IF ( .NOT. pmc_is_rootmodel() )  THEN
     742
    694743       CALL pmc_clientinit
    695744!
     
    725774
    726775       IF ( myid == 0 )  THEN
     776
    727777          CALL pmc_send_to_server( val, SIZE( val ), 0, 123, ierr )
    728778          CALL pmc_send_to_server( fval, SIZE( fval ), 0, 124, ierr )
     
    732782!
    733783!--       Receive Coarse grid information.
     784!--       TO_DO: find shorter and more meaningful name for  define_coarse_grid_real
    734785          CALL pmc_recv_from_server( define_coarse_grid_real,                  &
    735786                                     SIZE(define_coarse_grid_real), 0, 21, ierr )
     
    738789!
    739790!--       Receive also the dz-,zu- and zw-arrays here.
    740 !--       TO_DO: what is the meaning of above comment + remove write statements
    741 !--              and give this informations in header
    742           WRITE(0,*) 'Coarse grid from Server '
    743           WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
    744           WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
    745           WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(5)
    746           WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(6)
    747           WRITE(0,*) 'dx            = ',define_coarse_grid_real(3)
    748           WRITE(0,*) 'dy            = ',define_coarse_grid_real(4)
    749           WRITE(0,*) 'dz            = ',define_coarse_grid_real(7)
    750           WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
    751           WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
    752           WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)       
     791!--       TO_DO: what is the meaning of above comment
     792!
     793!--        Debug-printouts - keep them
     794!          WRITE(0,*) 'Coarse grid from Server '
     795!          WRITE(0,*) 'startx_tot    = ',define_coarse_grid_real(1)
     796!          WRITE(0,*) 'starty_tot    = ',define_coarse_grid_real(2)
     797!          WRITE(0,*) 'endx_tot      = ',define_coarse_grid_real(5)
     798!          WRITE(0,*) 'endy_tot      = ',define_coarse_grid_real(6)
     799!          WRITE(0,*) 'dx            = ',define_coarse_grid_real(3)
     800!          WRITE(0,*) 'dy            = ',define_coarse_grid_real(4)
     801!          WRITE(0,*) 'dz            = ',define_coarse_grid_real(7)
     802!          WRITE(0,*) 'nx_coarse     = ',define_coarse_grid_int(1)
     803!          WRITE(0,*) 'ny_coarse     = ',define_coarse_grid_int(2)
     804!          WRITE(0,*) 'nz_coarse     = ',define_coarse_grid_int(3)
    753805       ENDIF
    754806
     
    765817
    766818!
    767 !--    Get Server coordinates on coarse grid
     819!--    Get server coordinates on coarse grid
    768820       ALLOCATE( cg%coord_x(-nbgp:cg%nx+nbgp) )
    769821       ALLOCATE( cg%coord_y(-nbgp:cg%ny+nbgp) )
     
    776828!
    777829!--    Get coarse grid coordinates and vales of the z-direction from server
    778        IF ( myid == 0) THEN
    779           CALL pmc_recv_from_server( cg%coord_x, cg%nx + 1 + 2 * nbgp, 0, 24,  &
    780                                      ierr )
    781           CALL pmc_recv_from_server( cg%coord_y, cg%ny + 1 + 2 * nbgp, 0, 25,  &
    782                                      ierr )
     830       IF ( myid == 0)  THEN
     831
     832          CALL pmc_recv_from_server( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr )
     833          CALL pmc_recv_from_server( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr )
    783834          CALL pmc_recv_from_server( cg%dzu, cg%nz + 1, 0, 26, ierr )
    784835          CALL pmc_recv_from_server( cg%dzw, cg%nz + 1, 0, 27, ierr )
    785836          CALL pmc_recv_from_server( cg%zu, cg%nz + 2, 0, 28, ierr )
    786837          CALL pmc_recv_from_server( cg%zw, cg%nz + 2, 0, 29, ierr )
    787        ENDIF
    788 
    789 !
    790 !--    and broadcast this information
    791        CALL MPI_BCAST( cg%coord_x, cg%nx + 1 + 2 * nbgp, MPI_REAL, 0, comm2d,  &
    792                        ierr )
    793        CALL MPI_BCAST( cg%coord_y, cg%ny + 1 + 2 * nbgp, MPI_REAL, 0, comm2d,  &
    794                        ierr )
    795        CALL MPI_BCAST( cg%dzu, cg%nz + 1, MPI_REAL, 0, comm2d, ierr )
    796        CALL MPI_BCAST( cg%dzw, cg%nz + 1, MPI_REAL, 0, comm2d, ierr )
    797        CALL MPI_BCAST( cg%zu, cg%nz + 2,  MPI_REAL, 0, comm2d, ierr )
    798        CALL MPI_BCAST( cg%zw, cg%nz + 2,  MPI_REAL, 0, comm2d, ierr )
     838
     839       ENDIF
     840
     841!
     842!--    Broadcast this information
     843       CALL MPI_BCAST( cg%coord_x, cg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
     844       CALL MPI_BCAST( cg%coord_y, cg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
     845       CALL MPI_BCAST( cg%dzu, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
     846       CALL MPI_BCAST( cg%dzw, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
     847       CALL MPI_BCAST( cg%zu, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
     848       CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
    799849       
     850!
     851!--    TO_DO: give comments what is happening here
    800852       CALL pmci_map_fine_to_coarse_grid
    801853       CALL pmc_c_get_2d_index_list
    802854
    803855!
    804 !--    Include couple arrays into client content.
     856!--    Include couple arrays into client content
     857!--    TO_DO: Klaus: better explain the above comment (what is client content?)
    805858       CALL  pmc_c_clear_next_array_list
    806859       DO  WHILE ( pmc_c_getnextarray( myname ) )
     
    823876
    824877!
    825 !--    Two-way coupling     
     878!--    Two-way coupling
     879!--    TO_DO: comment what is happening here
    826880       IF ( nesting_mode == 'two-way' )  THEN
    827881          CALL pmci_init_anterp_tophat
     
    837891 CONTAINS
    838892
    839 
    840893    SUBROUTINE pmci_map_fine_to_coarse_grid
    841894
    842         IMPLICIT NONE
    843 
    844         INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all   !:
    845         INTEGER(iwp), DIMENSION(2)           ::  size_of_array      !:
     895       IMPLICIT NONE
     896
     897       INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all   !:
     898       INTEGER(iwp), DIMENSION(2)          ::  size_of_array      !:
    846899                                             
    847         REAL(wp)                             ::  coarse_dx   !:
    848         REAL(wp)                             ::  coarse_dy   !:
    849         REAL(wp)                             ::  loffset     !:
    850         REAL(wp)                             ::  noffset     !:
    851         REAL(wp)                             ::  roffset     !:
    852         REAL(wp)                             ::  soffset     !:
    853 
    854 !
    855 !--     Determine indices of interpolation/anterpolation area in coarse grid.
    856         coarse_dx = cg%dx
    857         coarse_dy = cg%dy
    858        
    859         loffset = MOD( coord_x(nxl), coarse_dx )       !  If the fine- and coarse grid nodes do not match.
    860         xexl    = coarse_dx + loffset
    861         nhll    = CEILING( xexl / coarse_dx )          !  This is needed in the anterpolation phase.
    862         xcs     = coord_x(nxl) - xexl
    863         DO  i = 0, cg%nx
    864            IF ( cg%coord_x(i) > xcs )  THEN
    865               icl = MAX( -1, i-1 ) 
    866               EXIT
    867            ENDIF
    868         ENDDO
    869 
    870         roffset = MOD( coord_x(nxr + 1), coarse_dx )   !  If the fine- and coarse grid nodes do not match.
    871         xexr    = coarse_dx + roffset
    872         nhlr    = CEILING( xexr / coarse_dx )          !  This is needed in the anterpolation phase.
    873         xce     = coord_x(nxr) + xexr
    874         DO  i = cg%nx, 0 , -1
    875            IF ( cg%coord_x(i) < xce )  THEN
    876               icr = MIN( cg%nx + 1, i + 1 )
    877               EXIT
    878            ENDIF
    879         ENDDO
    880 
    881         soffset = MOD( coord_y(nys), coarse_dy )       !  If the fine- and coarse grid nodes do not match
    882         yexs    = coarse_dy + soffset
    883         nhls    = CEILING( yexs / coarse_dy )          !  This is needed in the anterpolation phase.
    884         ycs     = coord_y(nys) - yexs
    885         DO  j = 0, cg%ny
    886            IF ( cg%coord_y(j) > ycs )  THEN
    887               jcs = MAX( -nbgp, j - 1 ) 
    888               EXIT
    889            ENDIF
    890         ENDDO
    891 
    892         noffset = MOD( coord_y(nyn + 1), coarse_dy )   !  If the fine- and coarse grid nodes do not match
    893         yexn    = coarse_dy + noffset
    894         nhln    = CEILING( yexn / coarse_dy )          !  This is needed in the anterpolation phase.
    895         yce     = coord_y(nyn) + yexn
    896         DO  j = cg%ny, 0, -1
    897            IF ( cg%coord_y(j) < yce )  THEN
    898               jcn = MIN( cg%ny + nbgp, j + 1 ) 
    899               EXIT
    900            ENDIF
    901         ENDDO
    902 
    903         coarse_bound(1) = icl
    904         coarse_bound(2) = icr
    905         coarse_bound(3) = jcs
    906         coarse_bound(4) = jcn
    907         coarse_bound(5) = myid
    908 !
    909 !--     Note that MPI_Gather receives data from all processes in the rank order.
    910         CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5, &
    911                          MPI_INTEGER, 0, comm2d, ierr )
    912 
    913         IF ( myid == 0 )  THEN
    914            size_of_array(1) = SIZE( coarse_bound_all, 1 )
    915            size_of_array(2) = SIZE( coarse_bound_all, 2 )
    916            CALL pmc_send_to_server( size_of_array, 2, 0, 40, ierr )
    917            CALL pmc_send_to_server( coarse_bound_all, SIZE( coarse_bound_all ), 0, 41, ierr )
    918         ENDIF
    919 
    920    END SUBROUTINE pmci_map_fine_to_coarse_grid
    921 
    922 
    923 
    924    SUBROUTINE pmci_init_interp_tril
    925 
    926 !
    927 !--   Precomputation of the interpolation coefficients and client-array indices
    928 !--   to be used by the interpolation routines interp_tril_lr, interp_tril_ns and
    929 !--   interp_tril_t. Constant dz is still assumed.
    930 !
    931 !--                                  Antti Hellsten 3.3.2015.
    932 !
    933 !--   Modified for variable dz, but not yet tested.
    934 !--                                  Antti Hellsten 27.3.2015.
    935 !
    936       IMPLICIT NONE
    937 
    938       INTEGER(iwp) ::  i       !:
    939       INTEGER(iwp) ::  i1      !:
    940       INTEGER(iwp) ::  j       !:
    941       INTEGER(iwp) ::  j1      !:
    942       INTEGER(iwp) ::  k       !:
    943       INTEGER(iwp) ::  kc      !:
    944 
    945       REAL(wp) ::  coarse_dx   !:
    946       REAL(wp) ::  coarse_dy   !:
    947       REAL(wp) ::  coarse_dz   !:
    948       REAL(wp) ::  xb          !:
    949       REAL(wp) ::  xcsu        !:
    950       REAL(wp) ::  xfso        !:
    951       REAL(wp) ::  xcso        !:
    952       REAL(wp) ::  xfsu        !:
    953       REAL(wp) ::  yb          !:
    954       REAL(wp) ::  ycso        !:
    955       REAL(wp) ::  ycsv        !:
    956       REAL(wp) ::  yfso        !:
    957       REAL(wp) ::  yfsv        !:
    958       REAL(wp) ::  zcso        !:
    959       REAL(wp) ::  zcsw        !:
    960       REAL(wp) ::  zfso        !:
    961       REAL(wp) ::  zfsw        !:
     900       REAL(wp) ::  loffset     !:
     901       REAL(wp) ::  noffset     !:
     902       REAL(wp) ::  roffset     !:
     903       REAL(wp) ::  soffset     !:
     904
     905!
     906!--    Determine indices of interpolation/anterpolation area in the coarse grid
     907!--    If the fine- and coarse grid nodes do not match.
     908       loffset = MOD( coord_x(nxl), cg%dx )
     909       xexl    = cg%dx + loffset
     910!
     911!--    This is needed in the anterpolation phase
     912       nhll = CEILING( xexl / cg%dx )
     913       xcs  = coord_x(nxl) - xexl
     914       DO  i = 0, cg%nx
     915          IF ( cg%coord_x(i) > xcs )  THEN
     916             icl = MAX( -1, i-1 )
     917             EXIT
     918          ENDIF
     919       ENDDO
     920!
     921!--    If the fine- and coarse grid nodes do not match
     922       roffset = MOD( coord_x(nxr+1), cg%dx )
     923       xexr    = cg%dx + roffset
     924!
     925!--    This is needed in the anterpolation phase
     926       nhlr = CEILING( xexr / cg%dx )
     927       xce  = coord_x(nxr) + xexr
     928       DO  i = cg%nx, 0 , -1
     929          IF ( cg%coord_x(i) < xce )  THEN
     930             icr = MIN( cg%nx+1, i+1 )
     931             EXIT
     932          ENDIF
     933       ENDDO
     934!
     935!--    If the fine- and coarse grid nodes do not match
     936       soffset = MOD( coord_y(nys), cg%dy )
     937       yexs    = cg%dy + soffset
     938!
     939!--    This is needed in the anterpolation phase
     940       nhls = CEILING( yexs / cg%dy )
     941       ycs  = coord_y(nys) - yexs
     942       DO  j = 0, cg%ny
     943          IF ( cg%coord_y(j) > ycs )  THEN
     944             jcs = MAX( -nbgp, j-1 )
     945             EXIT
     946          ENDIF
     947       ENDDO
     948!
     949!--    If the fine- and coarse grid nodes do not match
     950       noffset = MOD( coord_y(nyn+1), cg%dy )
     951       yexn    = cg%dy + noffset
     952!
     953!--    This is needed in the anterpolation phase
     954       nhln = CEILING( yexn / cg%dy )
     955       yce  = coord_y(nyn) + yexn
     956       DO  j = cg%ny, 0, -1
     957          IF ( cg%coord_y(j) < yce )  THEN
     958             jcn = MIN( cg%ny + nbgp, j+1 )
     959             EXIT
     960          ENDIF
     961       ENDDO
     962
     963       coarse_bound(1) = icl
     964       coarse_bound(2) = icr
     965       coarse_bound(3) = jcs
     966       coarse_bound(4) = jcn
     967       coarse_bound(5) = myid
     968!
     969!--    Note that MPI_Gather receives data from all processes in the rank order
     970!--    TO_DO: refer to the line where this fact becomes important
     971       CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5, &
     972                        MPI_INTEGER, 0, comm2d, ierr )
     973
     974       IF ( myid == 0 )  THEN
     975          size_of_array(1) = SIZE( coarse_bound_all, 1 )
     976          size_of_array(2) = SIZE( coarse_bound_all, 2 )
     977          CALL pmc_send_to_server( size_of_array, 2, 0, 40, ierr )
     978          CALL pmc_send_to_server( coarse_bound_all, SIZE( coarse_bound_all ), &
     979                                   0, 41, ierr )
     980       ENDIF
     981
     982    END SUBROUTINE pmci_map_fine_to_coarse_grid
     983
     984
     985
     986    SUBROUTINE pmci_init_interp_tril
     987!
     988!--    Precomputation of the interpolation coefficients and client-array indices
     989!--    to be used by the interpolation routines interp_tril_lr, interp_tril_ns
     990!--    and interp_tril_t. Constant dz is still assumed.
     991
     992       IMPLICIT NONE
     993
     994       INTEGER(iwp) ::  i       !:
     995       INTEGER(iwp) ::  i1      !:
     996       INTEGER(iwp) ::  j       !:
     997       INTEGER(iwp) ::  j1      !:
     998       INTEGER(iwp) ::  k       !:
     999       INTEGER(iwp) ::  kc      !:
     1000
     1001       REAL(wp) ::  xb          !:
     1002       REAL(wp) ::  xcsu        !:
     1003       REAL(wp) ::  xfso        !:
     1004       REAL(wp) ::  xcso        !:
     1005       REAL(wp) ::  xfsu        !:
     1006       REAL(wp) ::  yb          !:
     1007       REAL(wp) ::  ycso        !:
     1008       REAL(wp) ::  ycsv        !:
     1009       REAL(wp) ::  yfso        !:
     1010       REAL(wp) ::  yfsv        !:
     1011       REAL(wp) ::  zcso        !:
     1012       REAL(wp) ::  zcsw        !:
     1013       REAL(wp) ::  zfso        !:
     1014       REAL(wp) ::  zfsw        !:
    9621015     
    9631016
    964       coarse_dx = cg%dx
    965       coarse_dy = cg%dy
    966       coarse_dz = cg%dz
    967       xb        = nxl * dx
    968       yb        = nys * dy
     1017       xb = nxl * dx
     1018       yb = nys * dy
    9691019     
    970       ALLOCATE( icu(nxlg:nxrg) )
    971       ALLOCATE( ico(nxlg:nxrg) )
    972       ALLOCATE( jcv(nysg:nyng) )
    973       ALLOCATE( jco(nysg:nyng) )
    974       ALLOCATE( kcw(nzb:nzt+1) )
    975       ALLOCATE( kco(nzb:nzt+1) )
    976       ALLOCATE( r1xu(nxlg:nxrg) )
    977       ALLOCATE( r2xu(nxlg:nxrg) )
    978       ALLOCATE( r1xo(nxlg:nxrg) )
    979       ALLOCATE( r2xo(nxlg:nxrg) )
    980       ALLOCATE( r1yv(nysg:nyng) )
    981       ALLOCATE( r2yv(nysg:nyng) )
    982       ALLOCATE( r1yo(nysg:nyng) )
    983       ALLOCATE( r2yo(nysg:nyng) )
    984       ALLOCATE( r1zw(nzb:nzt+1) )
    985       ALLOCATE( r2zw(nzb:nzt+1) )
    986       ALLOCATE( r1zo(nzb:nzt+1) )
    987       ALLOCATE( r2zo(nzb:nzt+1) )
    988 
    989 !
    990 !--   Note that the node coordinates xfs... and xcs... are relative
    991 !--   to the lower-left-bottom corner of the fc-array, not the actual
    992 !--   client domain corner.
    993       DO  i = nxlg, nxrg
    994          xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
    995          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
    996          icu(i)  = icl + FLOOR( xfsu / coarse_dx )
    997          ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * coarse_dx ) / coarse_dx )
    998          xcsu    = ( icu(i) - icl ) * coarse_dx
    999          xcso    = ( ico(i) - icl ) * coarse_dx + 0.5_wp * coarse_dx
    1000          r2xu(i) = ( xfsu - xcsu ) / coarse_dx
    1001          r2xo(i) = ( xfso - xcso ) / coarse_dx
    1002          r1xu(i) = 1.0_wp - r2xu(i)
    1003          r1xo(i) = 1.0_wp - r2xo(i)
    1004       ENDDO
    1005 
    1006       DO  j = nysg, nyng
    1007          yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
    1008          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
    1009          jcv(j)  = jcs + FLOOR( yfsv/coarse_dy )           
    1010          jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * coarse_dy ) / coarse_dy )
    1011          ycsv    = ( jcv(j) - jcs ) * coarse_dy
    1012          ycso    = ( jco(j) - jcs ) * coarse_dy + 0.5_wp * coarse_dy         
    1013          r2yv(j) = ( yfsv - ycsv ) / coarse_dy           
    1014          r2yo(j) = ( yfso - ycso ) / coarse_dy 
    1015          r1yv(j) = 1.0_wp - r2yv(j)
    1016          r1yo(j) = 1.0_wp - r2yo(j)       
    1017       ENDDO
    1018 
    1019       DO  k = nzb, nzt + 1
    1020          zfsw = zw(k)
    1021          zfso = zu(k)
    1022 
    1023          kc = 0
    1024          DO  WHILE ( cg%zw(kc) <= zfsw )
    1025             kc = kc + 1
    1026          ENDDO
    1027          kcw(k) = kc - 1
     1020       ALLOCATE( icu(nxlg:nxrg) )
     1021       ALLOCATE( ico(nxlg:nxrg) )
     1022       ALLOCATE( jcv(nysg:nyng) )
     1023       ALLOCATE( jco(nysg:nyng) )
     1024       ALLOCATE( kcw(nzb:nzt+1) )
     1025       ALLOCATE( kco(nzb:nzt+1) )
     1026       ALLOCATE( r1xu(nxlg:nxrg) )
     1027       ALLOCATE( r2xu(nxlg:nxrg) )
     1028       ALLOCATE( r1xo(nxlg:nxrg) )
     1029       ALLOCATE( r2xo(nxlg:nxrg) )
     1030       ALLOCATE( r1yv(nysg:nyng) )
     1031       ALLOCATE( r2yv(nysg:nyng) )
     1032       ALLOCATE( r1yo(nysg:nyng) )
     1033       ALLOCATE( r2yo(nysg:nyng) )
     1034       ALLOCATE( r1zw(nzb:nzt+1) )
     1035       ALLOCATE( r2zw(nzb:nzt+1) )
     1036       ALLOCATE( r1zo(nzb:nzt+1) )
     1037       ALLOCATE( r2zo(nzb:nzt+1) )
     1038
     1039!
     1040!--    Note that the node coordinates xfs... and xcs... are relative to the
     1041!--    lower-left-bottom corner of the fc-array, not the actual client domain
     1042!--    corner
     1043       DO  i = nxlg, nxrg
     1044          xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
     1045          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
     1046          icu(i)  = icl + FLOOR( xfsu / cg%dx )
     1047          ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * cg%dx ) / cg%dx )
     1048          xcsu    = ( icu(i) - icl ) * cg%dx
     1049          xcso    = ( ico(i) - icl ) * cg%dx + 0.5_wp * cg%dx
     1050          r2xu(i) = ( xfsu - xcsu ) / cg%dx
     1051          r2xo(i) = ( xfso - xcso ) / cg%dx
     1052          r1xu(i) = 1.0_wp - r2xu(i)
     1053          r1xo(i) = 1.0_wp - r2xo(i)
     1054       ENDDO
     1055
     1056       DO  j = nysg, nyng
     1057          yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
     1058          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
     1059          jcv(j)  = jcs + FLOOR( yfsv / cg%dy )
     1060          jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * cg%dy ) / cg%dy )
     1061          ycsv    = ( jcv(j) - jcs ) * cg%dy
     1062          ycso    = ( jco(j) - jcs ) * cg%dy + 0.5_wp * cg%dy
     1063          r2yv(j) = ( yfsv - ycsv ) / cg%dy
     1064          r2yo(j) = ( yfso - ycso ) / cg%dy
     1065          r1yv(j) = 1.0_wp - r2yv(j)
     1066          r1yo(j) = 1.0_wp - r2yo(j)
     1067       ENDDO
     1068
     1069       DO  k = nzb, nzt + 1
     1070          zfsw = zw(k)
     1071          zfso = zu(k)
     1072
     1073          kc = 0
     1074          DO  WHILE ( cg%zw(kc) <= zfsw )
     1075             kc = kc + 1
     1076          ENDDO
     1077          kcw(k) = kc - 1
    10281078         
    1029          kc = 0
    1030          DO  WHILE ( cg%zu(kc) <= zfso )
    1031             kc = kc + 1
    1032          ENDDO
    1033          kco(k) = kc - 1
    1034 
    1035          zcsw    = cg%zw(kcw(k))
    1036          zcso    = cg%zu(kco(k))
    1037          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k) + 1 )
    1038          r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k) + 1 )
    1039          r1zw(k) = 1.0_wp - r2zw(k)
    1040          r1zo(k) = 1.0_wp - r2zo(k)
    1041       ENDDO
     1079          kc = 0
     1080          DO  WHILE ( cg%zu(kc) <= zfso )
     1081             kc = kc + 1
     1082          ENDDO
     1083          kco(k) = kc - 1
     1084
     1085          zcsw    = cg%zw(kcw(k))
     1086          zcso    = cg%zu(kco(k))
     1087          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kcw(k)+1)
     1088          r2zo(k) = ( zfso - zcso ) / cg%dzu(kco(k)+1)
     1089          r1zw(k) = 1.0_wp - r2zw(k)
     1090          r1zo(k) = 1.0_wp - r2zo(k)
     1091       ENDDO
    10421092     
    1043    END SUBROUTINE pmci_init_interp_tril
    1044 
    1045 
    1046 
    1047    SUBROUTINE pmci_init_loglaw_correction
    1048 
    1049 !
    1050 !--   Precomputation of the index and log-ratio arrays for the log-law corrections
    1051 !--   for near-wall nodes after the nest-BC interpolation.
    1052 !--   These are used by the interpolation routines interp_tril_lr and interp_tril_ns.
    1053 !
    1054 !--                                  Antti Hellsten 30.12.2015.
    1055 !
    1056       IMPLICIT NONE
    1057       INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i. 
    1058       INTEGER(iwp) ::  i            !:
    1059       INTEGER(iwp) ::  icorr        !:
    1060       INTEGER(iwp) ::  inc          !:  Wall outward-normal index increment -1 or 1, for direction=1, inc=1 always.
    1061       INTEGER(iwp) ::  iw           !:
    1062       INTEGER(iwp) ::  j            !:
    1063       INTEGER(iwp) ::  jcorr        !:
    1064       INTEGER(iwp) ::  jw           !:
    1065       INTEGER(iwp) ::  k            !:
    1066       INTEGER(iwp) ::  kb           !:
    1067       INTEGER(iwp) ::  kcorr        !:
    1068       INTEGER(iwp) ::  lc           !:
    1069       INTEGER(iwp) ::  ni           !:
    1070       INTEGER(iwp) ::  nj           !:
    1071       INTEGER(iwp) ::  nk           !:
    1072       INTEGER(iwp) ::  nzt_topo_max !:
    1073       INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
    1074 
    1075       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
    1076 
    1077 !
    1078 !--   First determine the maximum k-index needed for the near-wall corrections.
    1079 !--   This maximum is individual for each boundary to minimize the storage requirements
    1080 !--   and to minimize the corresponding loop k-range in the interpolation routines.
    1081       nzt_topo_nestbc_l = nzb
    1082       IF ( nest_bound_l )  THEN
    1083          DO  i = nxl - 1, nxl
    1084             DO  j = nys, nyn
    1085                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),   &
    1086                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
    1087             ENDDO
    1088          ENDDO
    1089          nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
    1090       ENDIF
     1093    END SUBROUTINE pmci_init_interp_tril
     1094
     1095
     1096
     1097    SUBROUTINE pmci_init_loglaw_correction
     1098!
     1099!--    Precomputation of the index and log-ratio arrays for the log-law
     1100!--    corrections for near-wall nodes after the nest-BC interpolation.
     1101!--    These are used by the interpolation routines interp_tril_lr and
     1102!--    interp_tril_ns.
     1103
     1104       IMPLICIT NONE
     1105
     1106       INTEGER(iwp) ::  direction    !:  Wall normal index: 1=k, 2=j, 3=i.
     1107       INTEGER(iwp) ::  i            !:
     1108       INTEGER(iwp) ::  icorr        !:
     1109       INTEGER(iwp) ::  inc          !:  Wall outward-normal index increment -1
     1110                                     !: or 1, for direction=1, inc=1 always
     1111       INTEGER(iwp) ::  iw           !:
     1112       INTEGER(iwp) ::  j            !:
     1113       INTEGER(iwp) ::  jcorr        !:
     1114       INTEGER(iwp) ::  jw           !:
     1115       INTEGER(iwp) ::  k            !:
     1116       INTEGER(iwp) ::  kb           !:
     1117       INTEGER(iwp) ::  kcorr        !:
     1118       INTEGER(iwp) ::  lc           !:
     1119       INTEGER(iwp) ::  ni           !:
     1120       INTEGER(iwp) ::  nj           !:
     1121       INTEGER(iwp) ::  nk           !:
     1122       INTEGER(iwp) ::  nzt_topo_max !:
     1123       INTEGER(iwp) ::  wall_index   !:  Index of the wall-node coordinate
     1124
     1125       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !:
     1126
     1127!
     1128!--    First determine the maximum k-index needed for the near-wall corrections.
     1129!--    This maximum is individual for each boundary to minimize the storage
     1130!--    requirements and to minimize the corresponding loop k-range in the
     1131!--    interpolation routines.
     1132       nzt_topo_nestbc_l = nzb
     1133       IF ( nest_bound_l )  THEN
     1134          DO  i = nxl-1, nxl
     1135             DO  j = nys, nyn
     1136                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l, nzb_u_inner(j,i),  &
     1137                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
     1138             ENDDO
     1139          ENDDO
     1140          nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
     1141       ENDIF
    10911142     
    1092       nzt_topo_nestbc_r = nzb
    1093       IF ( nest_bound_r )  THEN
    1094          i = nxr + 1
    1095          DO  j = nys, nyn
    1096             nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),      &
    1097                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
    1098          ENDDO
    1099          nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1       
    1100       ENDIF
    1101 
    1102       nzt_topo_nestbc_s = nzb
    1103       IF ( nest_bound_s )  THEN
    1104          DO  j = nys - 1, nys
    1105             DO  i = nxl, nxr
    1106                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),   &
    1107                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
    1108             ENDDO
    1109          ENDDO
    1110          nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1     
    1111       ENDIF
    1112 
    1113       nzt_topo_nestbc_n = nzb
    1114       IF ( nest_bound_n )  THEN
    1115          j = nyn + 1
    1116          DO  i = nxl, nxr
    1117             nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),      &
    1118                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
    1119          ENDDO
    1120          nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
    1121       ENDIF
    1122 
    1123 !
    1124 !--   Then determine the maximum number of near-wall nodes per wall point based
    1125 !--   on the grid-spacing ratios.
    1126       nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,                &
    1127                           nzt_topo_nestbc_s, nzt_topo_nestbc_n )
    1128 
    1129 !
    1130 !--   Note that the outer division must be integer division.
    1131       ni = CEILING( cg%dx / dx ) / 2
    1132       nj = CEILING( cg%dy / dy ) / 2
    1133       nk = 1
    1134       DO  k = 1, nzt_topo_max
    1135          nk = MAX( nk, CEILING( cg%dzu(k) / dzu(k) ) )
    1136       ENDDO
    1137       nk = nk / 2   !  Note that this must be integer division.
    1138       ncorr =  MAX( ni, nj, nk )
    1139 
    1140       ALLOCATE( lcr(0:ncorr - 1) )
    1141       lcr = 1.0_wp
    1142 
    1143 !
    1144 !--   First horizontal walls
    1145 !--   Left boundary
    1146       IF ( nest_bound_l )   THEN
    1147          ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
    1148          ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2) )
    1149          ALLOCATE( logc_ratio_u_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
    1150          ALLOCATE( logc_ratio_v_l(nzb:nzt_topo_nestbc_l, nys:nyn, 1:2, 0:ncorr-1) )
    1151          logc_u_l       = 0
    1152          logc_v_l       = 0
    1153          logc_ratio_u_l = 1.0_wp
    1154          logc_ratio_v_l = 1.0_wp
    1155          direction      = 1
    1156          inc            = 1
    1157 
    1158          DO  j = nys, nyn     
    1159 !
    1160 !--         Left boundary for u
    1161             i   = 0
    1162             kb  = nzb_u_inner(j,i)
    1163             k   = kb + 1
    1164             wall_index = kb
    1165             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1166                                      wall_index, z0(j,i), kb, direction, ncorr )
    1167             logc_u_l(k,j,1) = lc
    1168             logc_ratio_u_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
    1169             lcr(0:ncorr-1) = 1.0_wp
    1170 
    1171 !
    1172 !--         Left boundary for v
    1173             i   = -1
    1174             kb  =  nzb_v_inner(j,i)
    1175             k   =  kb + 1
    1176             wall_index = kb
    1177             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1178                                      wall_index, z0(j,i), kb, direction, ncorr )
    1179             logc_v_l(k,j,1) = lc
    1180             logc_ratio_v_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
    1181             lcr(0:ncorr-1) = 1.0_wp
    1182 
    1183          ENDDO
    1184       ENDIF
    1185 
    1186 !
    1187 !--   Right boundary
    1188       IF ( nest_bound_r )  THEN
    1189          ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
    1190          ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
    1191          ALLOCATE( logc_ratio_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
    1192          ALLOCATE( logc_ratio_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
    1193          logc_u_r       = 0
    1194          logc_v_r       = 0
    1195          logc_ratio_u_r = 1.0_wp
    1196          logc_ratio_v_r = 1.0_wp
    1197          direction      = 1
    1198          inc            = 1 
    1199          DO  j = nys, nyn
    1200 !
    1201 !--         Right boundary for u.
    1202             i   = nxr + 1
    1203             kb  = nzb_u_inner(j,i)
    1204             k   = kb + 1
    1205             wall_index = kb
    1206             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1207                                      wall_index, z0(j,i), kb, direction, ncorr )
    1208             logc_u_r(k,j,1) = lc
    1209             logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
    1210             lcr(0:ncorr-1) = 1.0_wp
    1211 
    1212 !
    1213 !--         Right boundary for v.
    1214             i   = nxr + 1
    1215             kb  = nzb_v_inner(j,i)
    1216             k   = kb + 1
    1217             wall_index = kb
    1218             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1219                                      wall_index, z0(j,i), kb, direction, ncorr )
    1220             logc_v_r(k,j,1) = lc
    1221             logc_ratio_v_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
    1222             lcr(0:ncorr-1) = 1.0_wp
    1223 
    1224          ENDDO
    1225       ENDIF
    1226 
    1227 !
    1228 !--   South boundary
    1229       IF ( nest_bound_s )  THEN
    1230          ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
    1231          ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
    1232          ALLOCATE( logc_ratio_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
    1233          ALLOCATE( logc_ratio_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
    1234          logc_u_s       = 0
    1235          logc_v_s       = 0
    1236          logc_ratio_u_s = 1.0_wp
    1237          logc_ratio_v_s = 1.0_wp
    1238          direction      = 1
    1239          inc            = 1
    1240          DO  i = nxl, nxr
    1241 !
    1242 !--         South boundary for u.
    1243             j   = -1
    1244             kb  =  nzb_u_inner(j,i)
    1245             k   =  kb + 1
    1246             wall_index = kb
    1247             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1248                                      wall_index, z0(j,i), kb, direction, ncorr )
    1249             logc_u_s(k,i,1) = lc
    1250             logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
    1251             lcr(0:ncorr-1) = 1.0_wp
    1252 
    1253 !
    1254 !--         South boundary for v
    1255             j   = 0
    1256             kb  = nzb_v_inner(j,i)
    1257             k   = kb + 1
    1258             wall_index = kb
    1259             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1260                                      wall_index, z0(j,i), kb, direction, ncorr )
    1261             logc_v_s(k,i,1) = lc
    1262             logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
    1263             lcr(0:ncorr-1) = 1.0_wp
    1264          ENDDO
    1265       ENDIF
    1266 
    1267 !
    1268 !--   North boundary
    1269       IF ( nest_bound_n )  THEN
    1270          ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
    1271          ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
    1272          ALLOCATE( logc_ratio_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
    1273          ALLOCATE( logc_ratio_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
    1274          logc_u_n       = 0
    1275          logc_v_n       = 0
    1276          logc_ratio_u_n = 1.0_wp
    1277          logc_ratio_v_n = 1.0_wp
    1278          direction      = 1
    1279          inc            = 1
    1280          DO  i = nxl, nxr
    1281 !
    1282 !--         North boundary for u.
    1283             j   = nyn + 1
    1284             kb  = nzb_u_inner(j,i)
    1285             k   = kb + 1
    1286             wall_index = kb
    1287             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1288                                      wall_index, z0(j,i), kb, direction, ncorr )
    1289             logc_u_n(k,i,1) = lc
    1290             logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
    1291             lcr(0:ncorr-1) = 1.0_wp
    1292 
    1293 !
    1294 !--         North boundary for v.
    1295             j   = nyn + 1
    1296             kb  = nzb_v_inner(j,i)
    1297             k   = kb + 1
    1298             wall_index = kb
    1299             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, &
    1300                                      wall_index, z0(j,i), kb, direction, ncorr )
    1301             logc_v_n(k,i,1) = lc
    1302             logc_ratio_v_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
    1303             lcr(0:ncorr-1) = 1.0_wp
    1304 
    1305          ENDDO
    1306       ENDIF
     1143       nzt_topo_nestbc_r = nzb
     1144       IF ( nest_bound_r )  THEN
     1145          i = nxr + 1
     1146          DO  j = nys, nyn
     1147             nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r, nzb_u_inner(j,i),     &
     1148                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
     1149          ENDDO
     1150          nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
     1151       ENDIF
     1152
     1153       nzt_topo_nestbc_s = nzb
     1154       IF ( nest_bound_s )  THEN
     1155          DO  j = nys-1, nys
     1156             DO  i = nxl, nxr
     1157                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s, nzb_u_inner(j,i),  &
     1158                                         nzb_v_inner(j,i), nzb_w_inner(j,i) )
     1159             ENDDO
     1160          ENDDO
     1161          nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1
     1162       ENDIF
     1163
     1164       nzt_topo_nestbc_n = nzb
     1165       IF ( nest_bound_n )  THEN
     1166          j = nyn + 1
     1167          DO  i = nxl, nxr
     1168             nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n, nzb_u_inner(j,i),     &
     1169                                      nzb_v_inner(j,i), nzb_w_inner(j,i) )
     1170          ENDDO
     1171          nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
     1172       ENDIF
     1173
     1174!
     1175!--    Then determine the maximum number of near-wall nodes per wall point based
     1176!--    on the grid-spacing ratios.
     1177       nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,               &
     1178                           nzt_topo_nestbc_s, nzt_topo_nestbc_n )
     1179
     1180!
     1181!--    Note that the outer division must be integer division.
     1182       ni = CEILING( cg%dx / dx ) / 2
     1183       nj = CEILING( cg%dy / dy ) / 2
     1184       nk = 1
     1185       DO  k = 1, nzt_topo_max
     1186          nk = MAX( nk, CEILING( cg%dzu(k) / dzu(k) ) )
     1187       ENDDO
     1188       nk = nk / 2   !  Note that this must be integer division.
     1189       ncorr =  MAX( ni, nj, nk )
     1190
     1191       ALLOCATE( lcr(0:ncorr-1) )
     1192       lcr = 1.0_wp
     1193
     1194!
     1195!--    First horizontal walls
     1196!--    Left boundary
     1197       IF ( nest_bound_l )  THEN
     1198          ALLOCATE( logc_u_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
     1199          ALLOCATE( logc_v_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
     1200          ALLOCATE( logc_ratio_u_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
     1201          ALLOCATE( logc_ratio_v_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
     1202          logc_u_l       = 0
     1203          logc_v_l       = 0
     1204          logc_ratio_u_l = 1.0_wp
     1205          logc_ratio_v_l = 1.0_wp
     1206          direction      = 1
     1207          inc            = 1
     1208
     1209          DO  j = nys, nyn
     1210!
     1211!--          Left boundary for u
     1212             i   = 0
     1213             kb  = nzb_u_inner(j,i)
     1214             k   = kb + 1
     1215             wall_index = kb
     1216             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1217                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1218             logc_u_l(k,j,1) = lc
     1219             logc_ratio_u_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
     1220             lcr(0:ncorr-1) = 1.0_wp
     1221!
     1222!--          Left boundary for v
     1223             i   = -1
     1224             kb  =  nzb_v_inner(j,i)
     1225             k   =  kb + 1
     1226             wall_index = kb
     1227             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1228                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1229             logc_v_l(k,j,1) = lc
     1230             logc_ratio_v_l(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
     1231             lcr(0:ncorr-1) = 1.0_wp
     1232
     1233          ENDDO
     1234       ENDIF
     1235
     1236!
     1237!--    Right boundary
     1238       IF ( nest_bound_r )  THEN
     1239          ALLOCATE( logc_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
     1240          ALLOCATE( logc_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2) )
     1241          ALLOCATE( logc_ratio_u_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
     1242          ALLOCATE( logc_ratio_v_r(nzb:nzt_topo_nestbc_r,nys:nyn,1:2,0:ncorr-1) )
     1243          logc_u_r       = 0
     1244          logc_v_r       = 0
     1245          logc_ratio_u_r = 1.0_wp
     1246          logc_ratio_v_r = 1.0_wp
     1247          direction      = 1
     1248          inc            = 1
     1249          DO  j = nys, nyn
     1250!
     1251!--          Right boundary for u.
     1252             i   = nxr + 1
     1253             kb  = nzb_u_inner(j,i)
     1254             k   = kb + 1
     1255             wall_index = kb
     1256             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1257                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1258             logc_u_r(k,j,1) = lc
     1259             logc_ratio_u_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
     1260             lcr(0:ncorr-1) = 1.0_wp
     1261
     1262!
     1263!--          Right boundary for v.
     1264             i   = nxr + 1
     1265             kb  = nzb_v_inner(j,i)
     1266             k   = kb + 1
     1267             wall_index = kb
     1268             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1269                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1270             logc_v_r(k,j,1) = lc
     1271             logc_ratio_v_r(k,j,1,0:ncorr-1) = lcr(0:ncorr-1)
     1272             lcr(0:ncorr-1) = 1.0_wp
     1273
     1274          ENDDO
     1275       ENDIF
     1276
     1277!
     1278!--    South boundary
     1279       IF ( nest_bound_s )  THEN
     1280          ALLOCATE( logc_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
     1281          ALLOCATE( logc_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2) )
     1282          ALLOCATE( logc_ratio_u_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
     1283          ALLOCATE( logc_ratio_v_s(nzb:nzt_topo_nestbc_s,nxl:nxr,1:2,0:ncorr-1) )
     1284          logc_u_s       = 0
     1285          logc_v_s       = 0
     1286          logc_ratio_u_s = 1.0_wp
     1287          logc_ratio_v_s = 1.0_wp
     1288          direction      = 1
     1289          inc            = 1
     1290          DO  i = nxl, nxr
     1291!
     1292!--          South boundary for u.
     1293             j   = -1
     1294             kb  =  nzb_u_inner(j,i)
     1295             k   =  kb + 1
     1296             wall_index = kb
     1297             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1298                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1299             logc_u_s(k,i,1) = lc
     1300             logc_ratio_u_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
     1301             lcr(0:ncorr-1) = 1.0_wp
     1302
     1303!
     1304!--          South boundary for v
     1305             j   = 0
     1306             kb  = nzb_v_inner(j,i)
     1307             k   = kb + 1
     1308             wall_index = kb
     1309             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1310                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1311             logc_v_s(k,i,1) = lc
     1312             logc_ratio_v_s(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
     1313             lcr(0:ncorr-1) = 1.0_wp
     1314          ENDDO
     1315       ENDIF
     1316
     1317!
     1318!--    North boundary
     1319       IF ( nest_bound_n )  THEN
     1320          ALLOCATE( logc_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
     1321          ALLOCATE( logc_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2) )
     1322          ALLOCATE( logc_ratio_u_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
     1323          ALLOCATE( logc_ratio_v_n(nzb:nzt_topo_nestbc_n,nxl:nxr,1:2,0:ncorr-1) )
     1324          logc_u_n       = 0
     1325          logc_v_n       = 0
     1326          logc_ratio_u_n = 1.0_wp
     1327          logc_ratio_v_n = 1.0_wp
     1328          direction      = 1
     1329          inc            = 1
     1330          DO  i = nxl, nxr
     1331!
     1332!--          North boundary for u.
     1333             j   = nyn + 1
     1334             kb  = nzb_u_inner(j,i)
     1335             k   = kb + 1
     1336             wall_index = kb
     1337             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1338                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1339             logc_u_n(k,i,1) = lc
     1340             logc_ratio_u_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
     1341             lcr(0:ncorr-1) = 1.0_wp
     1342
     1343!
     1344!--          North boundary for v.
     1345             j   = nyn + 1
     1346             kb  = nzb_v_inner(j,i)
     1347             k   = kb + 1
     1348             wall_index = kb
     1349             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j,     &
     1350                                inc, wall_index, z0(j,i), kb, direction, ncorr )
     1351             logc_v_n(k,i,1) = lc
     1352             logc_ratio_v_n(k,i,1,0:ncorr-1) = lcr(0:ncorr-1)
     1353             lcr(0:ncorr-1) = 1.0_wp
     1354
     1355          ENDDO
     1356       ENDIF
    13071357
    13081358!       
    1309 !--   Then vertical walls and corners if necessary.
    1310       IF ( topography /= 'flat' ) THEN       
    1311          kb = 0       ! kb is not used when direction > 1
     1359!--    Then vertical walls and corners if necessary.
     1360       IF ( topography /= 'flat' )  THEN
     1361          kb = 0       ! kb is not used when direction > 1
    13121362!       
    1313 !--      Left boundary
    1314          IF ( nest_bound_l )  THEN
    1315             ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
    1316             ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,0:ncorr-1) )
    1317             logc_w_l       = 0                   
    1318             logc_ratio_w_l = 1.0_wp
    1319             direction      = 2
    1320             DO  j = nys, nyn
    1321                DO  k = nzb, nzt_topo_nestbc_l
    1322 
    1323 !
    1324 !--               Wall for u on the south side, but not on the north side.
    1325                   i  = 0
    1326                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.     &
    1327                        ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )  THEN
    1328                      inc        =  1
    1329                      wall_index =  j
    1330                      CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
    1331 
    1332 !
    1333 !--                  The direction of the wall-normal index is stored as the sign of the logc-element.
    1334                      logc_u_l(k,j,2) = inc * lc
    1335                      logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
    1336                      lcr(0:ncorr-1) = 1.0_wp
    1337                   ENDIF
    1338 
    1339 !
    1340 !--               Wall for u on the north side, but not on the south side.
    1341                   i  = 0
     1363!--       Left boundary
     1364          IF ( nest_bound_l )  THEN
     1365             ALLOCATE( logc_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2) )
     1366             ALLOCATE( logc_ratio_w_l(nzb:nzt_topo_nestbc_l,nys:nyn,1:2,       &
     1367                                      0:ncorr-1) )
     1368             logc_w_l       = 0
     1369             logc_ratio_w_l = 1.0_wp
     1370             direction      = 2
     1371             DO  j = nys, nyn
     1372                DO  k = nzb, nzt_topo_nestbc_l
     1373
     1374!
     1375!--                Wall for u on the south side, but not on the north side
     1376                   i  = 0
     1377                   IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) .AND.        &
     1378                        ( nzb_u_outer(j,i) == nzb_u_outer(j-1,i) ) )           &
     1379                   THEN
     1380                      inc        =  1
     1381                      wall_index =  j
     1382                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
     1383                          k, j, inc, wall_index, z0(j,i), kb, direction, ncorr )
     1384!
     1385!--                   The direction of the wall-normal index is stored as the
     1386!--                   sign of the logc-element.
     1387                      logc_u_l(k,j,2) = inc * lc
     1388                      logc_ratio_u_l(k,j,2,0:ncorr-1) = lcr(0:ncorr-1)
     1389                      lcr(0:ncorr-1) = 1.0_wp
     1390                   ENDIF
     1391!
     1392!--                Wall for u on the north side, but not on the south side
     1393                   i  = 0
     1394!--                TO_DO: routine must be indentet by 1 space from here on,
     1395!--                       and long lines must be wrapped
    13421396                  IF ( ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) .AND.     &
    13431397                       ( nzb_u_outer(j,i) == nzb_u_outer(j+1,i) ) )  THEN
     
    16191673
    16201674
    1621    SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc, wall_index, z0_l, kb, direction, ncorr )
    1622       IMPLICIT NONE
    1623       INTEGER(iwp), INTENT(IN)  ::  direction                 !:
    1624       INTEGER(iwp), INTENT(IN)  ::  ij                        !:
    1625       INTEGER(iwp), INTENT(IN)  ::  inc                       !:
    1626       INTEGER(iwp), INTENT(IN)  ::  k                         !:
    1627       INTEGER(iwp), INTENT(IN)  ::  kb                        !:
    1628       INTEGER(iwp), INTENT(OUT) ::  lc                        !:
    1629       INTEGER(iwp), INTENT(IN)  ::  ncorr                     !:
    1630       INTEGER(iwp), INTENT(IN)  ::  wall_index                !:
    1631       REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !:
    1632       REAL(wp), INTENT(IN)      ::  z0_l                      !:
     1675    SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc,  &
     1676                                        wall_index, z0_l, kb, direction, ncorr )
     1677
     1678       IMPLICIT NONE
     1679
     1680       INTEGER(iwp), INTENT(IN)  ::  direction                 !:
     1681       INTEGER(iwp), INTENT(IN)  ::  ij                        !:
     1682       INTEGER(iwp), INTENT(IN)  ::  inc                       !:
     1683       INTEGER(iwp), INTENT(IN)  ::  k                         !:
     1684       INTEGER(iwp), INTENT(IN)  ::  kb                        !:
     1685       INTEGER(iwp), INTENT(OUT) ::  lc                        !:
     1686       INTEGER(iwp), INTENT(IN)  ::  ncorr                     !:
     1687       INTEGER(iwp), INTENT(IN)  ::  wall_index                !:
     1688
     1689       INTEGER(iwp) ::  alcorr       !:
     1690       INTEGER(iwp) ::  corr_index   !:
     1691       INTEGER(iwp) ::  lcorr        !:
     1692
     1693       LOGICAL      ::  more         !:
     1694
     1695       REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !:
     1696       REAL(wp), INTENT(IN)      ::  z0_l                      !:
    16331697     
    1634       INTEGER(iwp) ::  alcorr       !:
    1635       INTEGER(iwp) ::  corr_index   !:
    1636       INTEGER(iwp) ::  lcorr        !:
    1637       REAL(wp)     ::  logvelc1     !:
    1638       LOGICAL      ::  more         !:
     1698       REAL(wp)     ::  logvelc1     !:
    16391699     
    16401700
    1641       SELECT CASE ( direction )
    1642 
    1643          CASE (1)   !  k
    1644             more = .TRUE.
    1645             lcorr = 0
    1646             DO  WHILE ( more .AND. lcorr <= ncorr - 1 )
    1647                corr_index = k + lcorr
    1648                IF ( lcorr == 0 )  THEN
    1649                   CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
    1650                ENDIF
     1701       SELECT CASE ( direction )
     1702
     1703          CASE (1)   !  k
     1704             more = .TRUE.
     1705             lcorr = 0
     1706             DO  WHILE ( more .AND. lcorr <= ncorr-1 )
     1707                corr_index = k + lcorr
     1708                IF ( lcorr == 0 )  THEN
     1709                   CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
     1710                ENDIF
    16511711               
    1652                IF ( corr_index < lc )  THEN
    1653                   lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
    1654                   more = .TRUE.
    1655                ELSE
    1656                   lcr(lcorr) = 1.0
    1657                   more = .FALSE.
    1658                ENDIF
    1659                lcorr = lcorr + 1             
    1660             ENDDO
    1661 
    1662          CASE (2)   !  j
    1663             more = .TRUE.
    1664             lcorr  = 0
    1665             alcorr = 0
    1666             DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
    1667                corr_index = ij + lcorr   !  In this case (direction = 2) ij is j
    1668                IF ( lcorr == 0 )  THEN
    1669                   CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index, z0_l, inc )
    1670                ENDIF
    1671 
    1672 !
    1673 !--            The role of inc here is to make the comparison operation "<" valid in both directions.
    1674                IF ( inc * corr_index < inc * lc )  THEN
    1675                   lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy  &
    1676                               - coord_y(wall_index) ) / z0_l ) / logvelc1
    1677                   more = .TRUE.
    1678                ELSE
    1679                   lcr(alcorr) = 1.0_wp
    1680                   more = .FALSE.
    1681                ENDIF
    1682                lcorr  = lcorr + inc
    1683                alcorr = ABS( lcorr )
    1684             ENDDO
    1685 
    1686          CASE (3)   !  i         
    1687             more = .TRUE.
    1688             lcorr  = 0
    1689             alcorr = 0
    1690             DO  WHILE ( more  .AND.  alcorr <= ncorr - 1 )
    1691                corr_index = ij + lcorr   !  In this case (direction = 3) ij is i
    1692                IF ( lcorr == 0 )  THEN
    1693                   CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index, z0_l, inc )
    1694                ENDIF
    1695 
    1696 !
    1697 !--            The role of inc here is to make the comparison operation "<" valid in both directions.
    1698                IF ( inc * corr_index < inc * lc )  THEN
    1699                   lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx &
    1700                        - coord_x(wall_index) ) / z0_l ) / logvelc1
    1701                   more = .TRUE.
    1702                ELSE
    1703                   lcr(alcorr) = 1.0_wp
    1704                   more = .FALSE.
    1705                ENDIF
    1706                lcorr  = lcorr + inc
    1707                alcorr = ABS( lcorr )
    1708             ENDDO
    1709  
    1710       END SELECT
    1711 
    1712    END SUBROUTINE pmci_define_loglaw_correction_parameters
    1713 
    1714 
    1715 
    1716    SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
    1717 
    1718 !
    1719 !--   Finds the pivot node and te log-law factor for near-wall nodes for
    1720 !--   which the wall-parallel velocity components will be log-law corrected
    1721 !--   after interpolation. This subroutine is only for horizontal walls.
    1722 !
    1723 !--                               Antti Hellsten 30.12.2015
    1724       IMPLICIT NONE
    1725       REAL(wp),INTENT(OUT)  ::  logzc1     !:
    1726       REAL(wp), INTENT(IN)  ::  z0_l       !:
    1727       INTEGER(iwp), INTENT(IN)    ::  kb   !:
    1728       INTEGER(iwp), INTENT(OUT)   ::  lc   !:
     1712                IF ( corr_index < lc )  THEN
     1713                   lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
     1714                   more = .TRUE.
     1715                ELSE
     1716                   lcr(lcorr) = 1.0
     1717                   more = .FALSE.
     1718                ENDIF
     1719                lcorr = lcorr + 1
     1720             ENDDO
     1721
     1722          CASE (2)   !  j
     1723             more = .TRUE.
     1724             lcorr  = 0
     1725             alcorr = 0
     1726             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
     1727                corr_index = ij + lcorr   ! In this case (direction = 2) ij is j
     1728                IF ( lcorr == 0 )  THEN
     1729                   CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index,  &
     1730                                                z0_l, inc )
     1731                ENDIF
     1732
     1733!
     1734!--             The role of inc here is to make the comparison operation "<"
     1735!--             valid in both directions
     1736                IF ( inc * corr_index < inc * lc )  THEN
     1737                   lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy   &
     1738                                         - coord_y(wall_index) ) / z0_l )      &
     1739                                 / logvelc1
     1740                   more = .TRUE.
     1741                ELSE
     1742                   lcr(alcorr) = 1.0_wp
     1743                   more = .FALSE.
     1744                ENDIF
     1745                lcorr  = lcorr + inc
     1746                alcorr = ABS( lcorr )
     1747             ENDDO
     1748
     1749          CASE (3)   !  i
     1750             more = .TRUE.
     1751             lcorr  = 0
     1752             alcorr = 0
     1753             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
     1754                corr_index = ij + lcorr   ! In this case (direction = 3) ij is i
     1755                IF ( lcorr == 0 )  THEN
     1756                   CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index,  &
     1757                                                z0_l, inc )
     1758                ENDIF
     1759!
     1760!--             The role of inc here is to make the comparison operation "<"
     1761!--             valid in both directions
     1762                IF ( inc * corr_index < inc * lc )  THEN
     1763                   lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx   &
     1764                                         - coord_x(wall_index) ) / z0_l )      &
     1765                                 / logvelc1
     1766                   more = .TRUE.
     1767                ELSE
     1768                   lcr(alcorr) = 1.0_wp
     1769                   more = .FALSE.
     1770                ENDIF
     1771                lcorr  = lcorr + inc
     1772                alcorr = ABS( lcorr )
     1773             ENDDO
     1774
     1775       END SELECT
     1776
     1777    END SUBROUTINE pmci_define_loglaw_correction_parameters
     1778
     1779
     1780
     1781    SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
     1782!
     1783!--    Finds the pivot node and te log-law factor for near-wall nodes for
     1784!--    which the wall-parallel velocity components will be log-law corrected
     1785!--    after interpolation. This subroutine is only for horizontal walls.
     1786
     1787       IMPLICIT NONE
     1788
     1789       INTEGER(iwp), INTENT(IN)  ::  kb   !:
     1790       INTEGER(iwp), INTENT(OUT) ::  lc   !:
     1791
     1792       INTEGER(iwp) ::  kbc    !:
     1793       INTEGER(iwp) ::  k1     !:
     1794
     1795       REAL(wp),INTENT(OUT) ::  logzc1     !:
     1796       REAL(wp), INTENT(IN) ::  z0_l       !:
     1797
     1798       REAL(wp) ::  zuc1   !:
     1799
     1800
     1801       kbc = nzb + 1
     1802!
     1803!--    kbc is the first coarse-grid point above the surface
     1804       DO  WHILE ( cg%zu(kbc) < zu(kb) )
     1805          kbc = kbc + 1
     1806       ENDDO
     1807       zuc1  = cg%zu(kbc)
     1808       k1    = kb + 1
     1809       DO  WHILE ( zu(k1) < zuc1 )  !  Important: must be <, not <=
     1810          k1 = k1 + 1
     1811       ENDDO
     1812       logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
     1813       lc = k1
     1814
     1815    END SUBROUTINE pmci_find_logc_pivot_k
     1816
     1817
     1818
     1819    SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
     1820!
     1821!--    Finds the pivot node and te log-law factor for near-wall nodes for
     1822!--    which the wall-parallel velocity components will be log-law corrected
     1823!--    after interpolation. This subroutine is only for vertical walls on
     1824!--    south/north sides of the node.
     1825
     1826       IMPLICIT NONE
     1827
     1828       INTEGER(iwp), INTENT(IN)  ::  inc    !:  increment must be 1 or -1.
     1829       INTEGER(iwp), INTENT(IN)  ::  j      !:
     1830       INTEGER(iwp), INTENT(IN)  ::  jw     !:
     1831       INTEGER(iwp), INTENT(OUT) ::  lc     !:
     1832
     1833       INTEGER(iwp) ::  j1       !:
     1834
     1835       REAL(wp), INTENT(IN) ::  z0_l   !:
     1836
     1837       REAL(wp) ::  logyc1   !:
     1838       REAL(wp) ::  yc1      !:
     1839
     1840!
     1841!--    yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from
     1842!--    the wall
     1843       yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy
     1844!
     1845!--    j1 is the first fine-grid index further away from the wall than yc1
     1846       j1 = j
     1847!
     1848!--    Important: must be <, not <=
     1849       DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )
     1850          j1 = j1 + inc
     1851       ENDDO
     1852
     1853       logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
     1854       lc = j1
     1855
     1856    END SUBROUTINE pmci_find_logc_pivot_j
     1857
     1858
     1859
     1860    SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
     1861!
     1862!--    Finds the pivot node and the log-law factor for near-wall nodes for
     1863!--    which the wall-parallel velocity components will be log-law corrected
     1864!--    after interpolation. This subroutine is only for vertical walls on
     1865!--    south/north sides of the node.
     1866
     1867       IMPLICIT NONE
     1868
     1869       INTEGER(iwp), INTENT(IN)  ::  i      !:
     1870       INTEGER(iwp), INTENT(IN)  ::  inc    !: increment must be 1 or -1.
     1871       INTEGER(iwp), INTENT(IN)  ::  iw     !:
     1872       INTEGER(iwp), INTENT(OUT) ::  lc     !:
     1873
     1874       INTEGER(iwp) ::  i1       !:
     1875
     1876       REAL(wp), INTENT(IN) ::  z0_l   !:
     1877
     1878       REAL(wp) ::  logxc1   !:
     1879       REAL(wp) ::  xc1      !:
     1880
     1881!
     1882!--    xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from
     1883!--    the wall
     1884       xc1  = coord_x(iw) + 0.5_wp * inc * cg%dx
     1885!
     1886!--    i1 is the first fine-grid index futher away from the wall than xc1.
     1887       i1 = i
     1888!
     1889!--    Important: must be <, not <=
     1890       DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc * xc1 )
     1891          i1 = i1 + inc
     1892       ENDDO
    17291893     
    1730       REAL(wp)     ::  zuc1   !:
    1731       INTEGER(iwp) ::  kbc    !:
    1732       INTEGER(iwp) ::  k1     !:
    1733 
    1734 
    1735       kbc = nzb + 1
    1736       DO  WHILE ( cg%zu(kbc) < zu(kb) )   !  kbc is the first coarse-grid point above the surface.
    1737          kbc = kbc + 1
    1738       ENDDO
    1739       zuc1  = cg%zu(kbc)
    1740       k1    = kb + 1 
    1741       DO  WHILE ( zu(k1) < zuc1 )         !  Important: must be <, not <=
    1742          k1 = k1 + 1
    1743       ENDDO
    1744       logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
    1745       lc = k1
    1746 
    1747    END SUBROUTINE pmci_find_logc_pivot_k
    1748 
    1749 
    1750 
    1751    SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
    1752 
    1753 !
    1754 !--   Finds the pivot node and te log-law factor for near-wall nodes for
    1755 !--   which the wall-parallel velocity components will be log-law corrected
    1756 !--   after interpolation. This subroutine is only for vertical walls on
    1757 !--   south/north sides of the node.
    1758 !
    1759 !--                               Antti Hellsten 5.1.2016
    1760       IMPLICIT NONE
    1761       REAL(wp), INTENT(IN)      ::  z0_l   !:
    1762       INTEGER(iwp), INTENT(IN)  ::  inc    !:  increment must be 1 or -1.
    1763       INTEGER(iwp), INTENT(IN)  ::  j      !:
    1764       INTEGER(iwp), INTENT(IN)  ::  jw     !:
    1765       INTEGER(iwp), INTENT(OUT) ::  lc     !:
    1766      
    1767       REAL(wp)     ::  logyc1   !:
    1768       REAL(wp)     ::  yc1      !:
    1769       INTEGER(iwp) ::  j1       !:
    1770 
    1771 !
    1772 !--   yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from the wall.
    1773       yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy 
    1774 
    1775 !
    1776 !--   j1 is the first fine-grid index futher away from the wall than yc1.
    1777       j1 = j   
    1778       DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )   !  Important: must be <, not <=
    1779          j1 = j1 + inc
    1780       ENDDO
    1781 
    1782       logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
    1783       lc = j1
    1784 
    1785    END SUBROUTINE pmci_find_logc_pivot_j
    1786 
    1787 
    1788 
    1789    SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
    1790 
    1791 !
    1792 !--   Finds the pivot node and the log-law factor for near-wall nodes for
    1793 !--   which the wall-parallel velocity components will be log-law corrected
    1794 !--   after interpolation. This subroutine is only for vertical walls on
    1795 !--   south/north sides of the node.
    1796 !
    1797 !--                               Antti Hellsten 8.1.2016
    1798 
    1799       IMPLICIT NONE
    1800       REAL(wp), INTENT(IN)      ::  z0_l   !:
    1801       INTEGER(iwp), INTENT(IN)  ::  i      !:
    1802       INTEGER(iwp), INTENT(IN)  ::  inc    !: increment must be 1 or -1.
    1803       INTEGER(iwp), INTENT(IN)  ::  iw     !:
    1804       INTEGER(iwp), INTENT(OUT) ::  lc     !:
    1805 
    1806       REAL(wp)     ::  logxc1   !:
    1807       REAL(wp)     ::  xc1      !:
    1808       INTEGER(iwp) ::  i1       !:
    1809 
    1810 !
    1811 !--   xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from the wall.
    1812       xc1  = coord_x(iw) + 0.5_wp *inc * cg%dx 
    1813 
    1814 !
    1815 !--   i1 is the first fine-grid index futher away from the wall than xc1.
    1816       i1 = i   
    1817       DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc *xc1 )   !  Important: must be <, not <=
    1818          i1 = i1 + inc
    1819       ENDDO
    1820      
    1821       logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
    1822       lc = i1
    1823 
    1824    END SUBROUTINE pmci_find_logc_pivot_i
    1825 
    1826 
    1827 
     1894       logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
     1895       lc = i1
     1896
     1897    END SUBROUTINE pmci_find_logc_pivot_i
     1898
     1899
     1900!-- TO_DO:  indentation and wrap long lines from here on to the end of the file
    18281901   SUBROUTINE pmci_init_anterp_tophat
    18291902!
     
    18311904!--   corresponding coarse-grid array index and the
    18321905!--   Under-relaxation coefficients to be used by anterp_tophat.
    1833 !
    1834 !--                               Antti Hellsten 9.10.2015.
     1906
    18351907      IMPLICIT NONE
    1836       INTEGER(iwp) ::  i        !:
     1908
     1909      INTEGER(iwp) ::  i        !: Fine-grid index
     1910      INTEGER(iwp) ::  ii       !: Coarse-grid index
    18371911      INTEGER(iwp) ::  istart   !:
    1838       INTEGER(iwp) ::  j        !:
     1912      INTEGER(iwp) ::  j        !: Fine-grid index
     1913      INTEGER(iwp) ::  jj       !: Coarse-grid index
    18391914      INTEGER(iwp) ::  jstart   !:
    1840       INTEGER(iwp) ::  k        !:
     1915      INTEGER(iwp) ::  k        !: Fine-grid index
     1916      INTEGER(iwp) ::  kk       !: Coarse-grid index
    18411917      INTEGER(iwp) ::  kstart   !:
    1842       INTEGER(iwp) ::  l        !:
    1843       INTEGER(iwp) ::  m        !:
    1844       INTEGER(iwp) ::  n        !:
    18451918      REAL(wp)     ::  xi       !:
    18461919      REAL(wp)     ::  eta      !:
    18471920      REAL(wp)     ::  zeta     !:
    18481921     
     1922
    18491923!
    18501924!--   Default values:
     
    18661940
    18671941!
    1868 !--   First determine kceu and kcew that are the coarse-grid upper bounds for index k.
    1869       n = 0
    1870       DO  WHILE ( cg%zu(n) < zu(nzt) )
    1871          n = n + 1
     1942!--   First determine kctu and kctw that are the coarse-grid upper bounds for index k.
     1943      kk = 0
     1944      DO  WHILE ( cg%zu(kk) < zu(nzt) )
     1945         kk = kk + 1
    18721946      ENDDO
    1873       kceu = n - 1
    1874 
    1875       n = 0
    1876       DO  WHILE ( cg%zw(n) < zw(nzt-1) )
    1877          n = n + 1
     1947      kctu = kk - 1
     1948
     1949      kk = 0
     1950      DO  WHILE ( cg%zw(kk) < zw(nzt-1) )
     1951         kk = kk + 1
    18781952      ENDDO
    1879       kcew = n - 1
     1953      kctw = kk - 1
    18801954
    18811955      ALLOCATE( iflu(icl:icr) )
     
    18871961      ALLOCATE( jfuv(jcs:jcn) )
    18881962      ALLOCATE( jfuo(jcs:jcn) )
    1889       ALLOCATE( kflw(0:kcew) )
    1890       ALLOCATE( kflo(0:kceu) )
    1891       ALLOCATE( kfuw(0:kcew) )
    1892       ALLOCATE( kfuo(0:kceu) )
     1963      ALLOCATE( kflw(0:kctw) )
     1964      ALLOCATE( kflo(0:kctu) )
     1965      ALLOCATE( kfuw(0:kctw) )
     1966      ALLOCATE( kfuo(0:kctu) )
    18931967
    18941968!
    18951969!--   i-indices of u for each l-index value.   
    18961970      istart = nxlg
    1897       DO  l = icl, icr 
     1971      DO  ii = icl, icr 
    18981972         i = istart
    1899          DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) )
     1973         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) )
    19001974            i = i + 1
    19011975         ENDDO
    1902          iflu(l) = MIN( MAX( i, nxlg ), nxrg )
    1903          DO  WHILE ( ( coord_x(i)  <  cg%coord_x(l) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) )
     1976         iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
     1977         DO  WHILE ( ( coord_x(i)  <  cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.  ( i < nxrg ) )
    19041978            i = i + 1
    19051979         ENDDO
    1906          ifuu(l) = MIN( MAX( i, nxlg ), nxrg )
    1907          istart = iflu(l)
     1980         ifuu(ii) = MIN( MAX( i, nxlg ), nxrg )
     1981         istart = iflu(ii)
    19081982      ENDDO
    19091983
     
    19111985!--   i-indices of others for each l-index value.
    19121986      istart = nxlg
    1913       DO  l = icl, icr   
     1987      DO  ii = icl, icr   
    19141988         i = istart
    1915          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) )  .AND.  ( i < nxrg ) )
     1989         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(ii) )  .AND.  ( i < nxrg ) )
    19161990            i = i + 1
    19171991         ENDDO
    1918          iflo(l) = MIN( MAX( i, nxlg ), nxrg )
    1919          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(l) + cg%dx )  .AND.  ( i < nxrg ) )
     1992         iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
     1993         DO  WHILE ( ( coord_x(i) + 0.5_wp * dx  <  cg%coord_x(ii) + cg%dx )  .AND.  ( i < nxrg ) )
    19201994            i = i + 1
    19211995         ENDDO
    1922          ifuo(l) = MIN(MAX(i,nxlg),nxrg)
    1923          istart = iflo(l)
     1996         ifuo(ii) = MIN(MAX(i,nxlg),nxrg)
     1997         istart = iflo(ii)
    19241998      ENDDO
    19251999
     
    19272001!--   j-indices of v for each m-index value.
    19282002      jstart = nysg
    1929       DO  m = jcs, jcn
     2003      DO  jj = jcs, jcn
    19302004         j = jstart
    1931          DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) )
     2005         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  ( j < nyng ) )
    19322006            j = j + 1
    19332007         ENDDO
    1934          jflv(m) = MIN( MAX( j, nysg ), nyng )
    1935          DO  WHILE ( ( coord_y(j)  <  cg%coord_y(m) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) )
     2008         jflv(jj) = MIN( MAX( j, nysg ), nyng )
     2009         DO  WHILE ( ( coord_y(j)  <  cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.  ( j < nyng ) )
    19362010            j = j + 1
    19372011         ENDDO
    1938          jfuv(m) = MIN( MAX( j, nysg ), nyng )
    1939          jstart = jflv(m)
     2012         jfuv(jj) = MIN( MAX( j, nysg ), nyng )
     2013         jstart = jflv(jj)
    19402014      ENDDO
    19412015
     
    19432017!--   j-indices of others for each m-index value.
    19442018      jstart = nysg
    1945       DO  m = jcs, jcn
     2019      DO  jj = jcs, jcn
    19462020         j = jstart
    1947          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) )  .AND.  ( j < nyng ) ) 
     2021         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(jj) )  .AND.  ( j < nyng ) ) 
    19482022            j = j + 1
    19492023         ENDDO
    1950          jflo(m) = MIN( MAX( j, nysg ), nyng )
    1951          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(m) + cg%dy )  .AND.  ( j < nyng ) )
     2024         jflo(jj) = MIN( MAX( j, nysg ), nyng )
     2025         DO  WHILE ( ( coord_y(j) + 0.5_wp * dy  <  cg%coord_y(jj) + cg%dy )  .AND.  ( j < nyng ) )
    19522026            j = j + 1
    19532027         ENDDO
    1954          jfuo(m) = MIN( MAX( j, nysg ), nyng )
    1955          jstart = jflv(m)
     2028         jfuo(jj) = MIN( MAX( j, nysg ), nyng )
     2029         jstart = jflv(jj)
    19562030      ENDDO
    19572031
     
    19612035      kflw(0) = 0
    19622036      kfuw(0) = 0
    1963       DO  n = 1, kcew
     2037      DO  kk = 1, kctw
    19642038         k = kstart
    1965          DO  WHILE ( ( zw(k) < cg%zw(n) - 0.5_wp * cg%dzw(n) )  .AND.  ( k < nzt ) )
     2039         DO  WHILE ( ( zw(k) < cg%zw(kk) - 0.5_wp * cg%dzw(kk) )  .AND.  ( k < nzt ) )
    19662040            k = k + 1
    19672041         ENDDO
    1968          kflw(n) = MIN( MAX( k, 1 ), nzt + 1 )
    1969          DO  WHILE ( ( zw(k) < cg%zw(n) + 0.5_wp * cg%dzw(n+1) )  .AND.  ( k < nzt ) )
     2042         kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
     2043         DO  WHILE ( ( zw(k) < cg%zw(kk) + 0.5_wp * cg%dzw(kk+1) )  .AND.  ( k < nzt ) )
    19702044            k = k + 1
    19712045         ENDDO
    1972          kfuw(n) = MIN( MAX( k, 1 ), nzt + 1 )
    1973          kstart = kflw(n)
     2046         kfuw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
     2047         kstart = kflw(kk)
    19742048      ENDDO
    19752049
     
    19792053      kflo(0) = 0
    19802054      kfuo(0) = 0
    1981       DO  n = 1, kceu
     2055      DO  kk = 1, kctu
    19822056         k = kstart
    1983          DO  WHILE ( ( zu(k) < cg%zu(n) - 0.5_wp * cg%dzu(n) )  .AND.  ( k < nzt ) )
     2057         DO  WHILE ( ( zu(k) < cg%zu(kk) - 0.5_wp * cg%dzu(kk) )  .AND.  ( k < nzt ) )
    19842058            k = k + 1
    19852059         ENDDO
    1986          kflo(n) = MIN( MAX( k, 1 ), nzt + 1 )
    1987          DO  WHILE ( ( zu(k) < cg%zu(n) + 0.5_wp * cg%dzu(n+1) )  .AND.  ( k < nzt ) ) 
     2060         kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
     2061         DO  WHILE ( ( zu(k) < cg%zu(kk) + 0.5_wp * cg%dzu(kk+1) )  .AND.  ( k < nzt ) ) 
    19882062            k = k + 1
    19892063         ENDDO
    1990          kfuo(n) = MIN( MAX( k-1, 1 ), nzt + 1 )
    1991          kstart = kflo(n)
     2064         kfuo(kk) = MIN( MAX( k-1, 1 ), nzt + 1 )
     2065         kstart = kflo(kk)
    19922066      ENDDO
    19932067     
     
    19952069!--   Spatial under-relaxation coefficients
    19962070      ALLOCATE( frax(icl:icr) )
    1997       DO  l = icl, icr
     2071
     2072      DO  ii = icl, icr
    19982073         IF ( nest_bound_l ) THEN
    1999             xi   = ( ( cg%coord_x(l) - lower_left_coord_x ) / anterp_relax_length_l )**4
     2074            xi   = ( MAX( 0.0_wp, ( cg%coord_x(ii) - lower_left_coord_x ) ) / anterp_relax_length_l )**4
    20002075         ELSEIF ( nest_bound_r ) THEN
    2001             xi   = ( ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(l) ) / anterp_relax_length_r )**4
     2076            xi   = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx - cg%coord_x(ii) ) ) / anterp_relax_length_r )**4
    20022077         ELSE
    20032078            xi   = 999999.9_wp
    20042079         ENDIF
    2005          frax(l)  = xi / ( 1.0_wp + xi )
     2080         frax(ii)  = xi / ( 1.0_wp + xi )
    20062081      ENDDO
    20072082
    20082083      ALLOCATE( fray(jcs:jcn) )
    2009       DO  m = jcs, jcn
     2084
     2085      DO  jj = jcs, jcn
    20102086         IF ( nest_bound_s ) THEN           
    2011             eta  = ( ( cg%coord_y(m) - lower_left_coord_y ) / anterp_relax_length_s )**4
     2087            eta  = ( MAX( 0.0_wp, ( cg%coord_y(jj) - lower_left_coord_y ) ) / anterp_relax_length_s )**4
    20122088         ELSEIF ( nest_bound_n ) THEN
    2013             eta  = ( (lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(m)) / anterp_relax_length_n )**4
     2089            eta  = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy - cg%coord_y(jj)) ) / anterp_relax_length_n )**4
    20142090         ELSE
    20152091            eta  = 999999.9_wp
    20162092         ENDIF
    2017          fray(m)  = eta / ( 1.0_wp + eta )
     2093         fray(jj)  = eta / ( 1.0_wp + eta )
    20182094      ENDDO
    20192095     
    2020       ALLOCATE( fraz(0:kceu) )
    2021       DO  n = 0, kceu       
    2022          zeta = ( ( zu(nzt) - cg%zu(n) ) / anterp_relax_length_t )**4       
    2023          fraz(n) = zeta / ( 1.0_wp + zeta )
     2096      ALLOCATE( fraz(0:kctu) )
     2097      DO  kk = 0, kctu       
     2098         zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4       
     2099         fraz(kk) = zeta / ( 1.0_wp + zeta )
    20242100      ENDDO
    20252101
     
    21732249
    21742250
     2251
     2252 SUBROUTINE pmci_set_array_pointer( name, client_id, nz_cl )
     2253
     2254    IMPLICIT NONE
     2255
     2256    INTEGER, INTENT(IN)          ::  client_id   !:
     2257    INTEGER, INTENT(IN)          ::  nz_cl       !:
     2258    CHARACTER(LEN=*), INTENT(IN) ::  name        !:
     2259
     2260#if defined( __parallel )
     2261    INTEGER(iwp) ::  ierr        !:
     2262    INTEGER(iwp) ::  istat       !:
     2263
     2264    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d        !:
     2265    REAL(wp), POINTER, DIMENSION(:,:)   ::  p_2d_sec    !:
     2266    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d        !:
     2267    REAL(wp), POINTER, DIMENSION(:,:,:) ::  p_3d_sec    !:
     2268
     2269
     2270    NULLIFY( p_3d )
     2271    NULLIFY( p_2d )
     2272
     2273!
     2274!-- List of array names, which can be coupled.
     2275!-- In case of 3D please change also the second array for the pointer version
     2276    IF ( TRIM(name) == "u"  )  p_3d => u
     2277    IF ( TRIM(name) == "v"  )  p_3d => v
     2278    IF ( TRIM(name) == "w"  )  p_3d => w
     2279    IF ( TRIM(name) == "e"  )  p_3d => e
     2280    IF ( TRIM(name) == "pt" )  p_3d => pt
     2281    IF ( TRIM(name) == "q"  )  p_3d => q
     2282!
     2283!-- Next line is just an example for a 2D array (not active for coupling!)
     2284!-- Please note, that z0 has to be declared as TARGET array in modules.f90
     2285!    IF ( TRIM(name) == "z0" )    p_2d => z0
     2286
     2287#if defined( __nopointer )
     2288    IF ( ASSOCIATED( p_3d ) )  THEN
     2289       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz )
     2290    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
     2291       CALL pmc_s_set_dataarray( client_id, p_2d )
     2292    ELSE
     2293!
     2294!--    Give only one message for the root domain
     2295       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
     2296
     2297          message_string = 'pointer for array "' // TRIM( name ) //            &
     2298                           '" can''t be associated'
     2299          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
     2300       ELSE
     2301!
     2302!--       Avoid others to continue
     2303          CALL MPI_BARRIER( comm2d, ierr )
     2304       ENDIF
     2305    ENDIF
     2306#else
     2307    IF ( TRIM(name) == "u"  )  p_3d_sec => u_2
     2308    IF ( TRIM(name) == "v"  )  p_3d_sec => v_2
     2309    IF ( TRIM(name) == "w"  )  p_3d_sec => w_2
     2310    IF ( TRIM(name) == "e"  )  p_3d_sec => e_2
     2311    IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
     2312    IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
     2313
     2314    IF ( ASSOCIATED( p_3d ) )  THEN
     2315       CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz, &
     2316                                 array_2 = p_3d_sec )
     2317    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
     2318       CALL pmc_s_set_dataarray( client_id, p_2d )
     2319    ELSE
     2320!
     2321!--    Give only one message for the root domain
     2322       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
     2323
     2324          message_string = 'pointer for array "' // TRIM( name ) //            &
     2325                           '" can''t be associated'
     2326          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
     2327       ELSE
     2328!
     2329!--       Avoid others to continue
     2330          CALL MPI_BARRIER( comm2d, ierr )
     2331       ENDIF
     2332
     2333   ENDIF
     2334#endif
     2335
     2336#endif
     2337 END SUBROUTINE pmci_set_array_pointer
     2338
     2339
     2340
     2341 SUBROUTINE pmci_create_client_arrays( name, is, ie, js, je, nzc  )
     2342
     2343    IMPLICIT NONE
     2344
     2345    CHARACTER(LEN=*), INTENT(IN) ::  name    !:
     2346
     2347    INTEGER(iwp), INTENT(IN) ::  ie      !:
     2348    INTEGER(iwp), INTENT(IN) ::  is      !:
     2349    INTEGER(iwp), INTENT(IN) ::  je      !:
     2350    INTEGER(iwp), INTENT(IN) ::  js      !:
     2351    INTEGER(iwp), INTENT(IN) ::  nzc     !:  Note that nzc is cg%nz
     2352
     2353#if defined( __parallel )
     2354    INTEGER(iwp) ::  ierr    !:
     2355    INTEGER(iwp) ::  istat   !:
     2356
     2357    REAL(wp), POINTER,DIMENSION(:,:)   ::  p_2d    !:
     2358    REAL(wp), POINTER,DIMENSION(:,:,:) ::  p_3d    !:
     2359
     2360
     2361    NULLIFY( p_3d )
     2362    NULLIFY( p_2d )
     2363
     2364!
     2365!-- List of array names, which can be coupled.
     2366!-- TO_DO: Antti: what is the meaning of the next line?
     2367!-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1.
     2368    IF ( TRIM( name ) == "u" )  THEN
     2369       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
     2370       p_3d => uc
     2371    ELSEIF ( TRIM( name ) == "v" )  THEN
     2372       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1, js:je, is:ie) )
     2373       p_3d => vc
     2374    ELSEIF ( TRIM( name ) == "w" )  THEN
     2375       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1, js:je, is:ie) )
     2376       p_3d => wc
     2377    ELSEIF ( TRIM( name ) == "e" )  THEN
     2378       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1, js:je, is:ie) )
     2379       p_3d => ec
     2380    ELSEIF ( TRIM( name ) == "pt")  THEN
     2381       IF ( .NOT. ALLOCATED( ptc ) ) ALLOCATE( ptc(0:nzc+1, js:je, is:ie) )
     2382       p_3d => ptc
     2383    ELSEIF ( TRIM( name ) == "q")  THEN
     2384       IF ( .NOT. ALLOCATED( qc ) ) ALLOCATE( qc(0:nzc+1, js:je, is:ie) )
     2385       p_3d => qc
     2386    !ELSEIF (trim(name) == "z0") then
     2387       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
     2388       !p_2d => z0c
     2389    ENDIF
     2390
     2391    IF ( ASSOCIATED( p_3d ) )  THEN
     2392       CALL pmc_c_set_dataarray( p_3d )
     2393    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
     2394       CALL pmc_c_set_dataarray( p_2d )
     2395    ELSE
     2396!
     2397!--    Give only one message for the first client domain
     2398       IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
     2399
     2400          message_string = 'pointer for array "' // TRIM( name ) //            &
     2401                           '" can''t be associated'
     2402          CALL message( 'pmci_create_client_arrays', 'PA0170', 3, 2, 0, 6, 0 )
     2403       ELSE
     2404!
     2405!--       Avoid others to continue
     2406          CALL MPI_BARRIER( comm2d, ierr )
     2407       ENDIF
     2408    ENDIF
     2409
     2410#endif
     2411 END SUBROUTINE pmci_create_client_arrays
     2412
     2413
     2414
     2415 SUBROUTINE pmci_server_initialize
     2416!-- TO_DO: add general explanations about what this subroutine does
     2417#if defined( __parallel )
     2418    IMPLICIT NONE
     2419
     2420    INTEGER(iwp) ::  client_id   !:
     2421    INTEGER(iwp) ::  m           !:
     2422
     2423    REAL(wp) ::  waittime    !:
     2424
     2425
     2426    DO  m = 1, SIZE( pmc_server_for_client ) - 1
     2427       client_id = pmc_server_for_client(m)
     2428       CALL pmc_s_fillbuffer( client_id, waittime=waittime )
     2429    ENDDO
     2430
     2431#endif
     2432 END SUBROUTINE pmci_server_initialize
     2433
     2434
     2435
     2436 SUBROUTINE pmci_client_initialize
     2437!-- TO_DO: add general explanations about what this subroutine does
     2438#if defined( __parallel )
     2439    IMPLICIT NONE
     2440
     2441    INTEGER(iwp) ::  i          !:
     2442    INTEGER(iwp) ::  icl        !:
     2443    INTEGER(iwp) ::  icr        !:
     2444    INTEGER(iwp) ::  j          !:
     2445    INTEGER(iwp) ::  jcn        !:
     2446    INTEGER(iwp) ::  jcs        !:
     2447
     2448    REAL(wp) ::  waittime   !:
     2449
     2450!
     2451!-- Root id is never a client
     2452    IF ( cpl_id > 1 )  THEN
     2453
     2454!
     2455!--    Client domain boundaries in the server index space
     2456       icl = coarse_bound(1)
     2457       icr = coarse_bound(2)
     2458       jcs = coarse_bound(3)
     2459       jcn = coarse_bound(4)
     2460
     2461!
     2462!--    Get data from server
     2463       CALL pmc_c_getbuffer( waittime = waittime )
     2464
     2465!
     2466!--    The interpolation.
     2467       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
     2468                                   r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
     2469       CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,   &
     2470                                   r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
     2471       CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,   &
     2472                                   r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
     2473       CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
     2474                                   r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
     2475       CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
     2476                                   r2yo, r1zo, r2zo, nzb_s_inner, 's' )
     2477       IF ( humidity  .OR.  passive_scalar )  THEN
     2478          CALL pmci_interp_tril_all ( q, qc, ico, jco, kco, r1xo, r2xo, r1yo,  &
     2479                                      r2yo, r1zo, r2zo, nzb_s_inner, 's' )
     2480       ENDIF
     2481
     2482       IF ( topography /= 'flat' )  THEN
     2483!
     2484!--       Inside buildings set velocities and TKE back to zero.
     2485!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
     2486!--       maybe revise later.
     2487          DO   i = nxlg, nxrg
     2488             DO   j = nysg, nyng
     2489                u(nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
     2490                v(nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
     2491                w(nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
     2492                e(nzb:nzb_s_inner(j,i),j,i)   = 0.0_wp
     2493                u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
     2494                v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
     2495                w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
     2496                e_p(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
     2497             ENDDO
     2498          ENDDO
     2499       ENDIF
     2500    ENDIF
     2501
     2502
     2503 CONTAINS
     2504
     2505
     2506    SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
     2507                                     r1z, r2z, kb, var )
     2508!
     2509!--    Interpolation of the internal values for the client-domain initialization
     2510!--    This subroutine is based on trilinear interpolation.
     2511!--    Coding based on interp_tril_lr/sn/t
     2512       IMPLICIT NONE
     2513
     2514       CHARACTER(LEN=1), INTENT(IN) :: var  !:
     2515
     2516       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !:
     2517       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !:
     2518       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !:
     2519       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !:
     2520
     2521       INTEGER(iwp) ::  i      !:
     2522       INTEGER(iwp) ::  ib     !:
     2523       INTEGER(iwp) ::  ie     !:
     2524       INTEGER(iwp) ::  j      !:
     2525       INTEGER(iwp) ::  jb     !:
     2526       INTEGER(iwp) ::  je     !:
     2527       INTEGER(iwp) ::  k      !:
     2528       INTEGER(iwp) ::  k1     !:
     2529       INTEGER(iwp) ::  kbc    !:
     2530       INTEGER(iwp) ::  l      !:
     2531       INTEGER(iwp) ::  m      !:
     2532       INTEGER(iwp) ::  n      !:
     2533
     2534       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
     2535       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc    !:
     2536       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x   !:
     2537       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x   !:
     2538       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y   !:
     2539       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y   !:
     2540       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z   !:
     2541       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z   !:
     2542
     2543       REAL(wp) ::  fk         !:
     2544       REAL(wp) ::  fkj        !:
     2545       REAL(wp) ::  fkjp       !:
     2546       REAL(wp) ::  fkp        !:
     2547       REAL(wp) ::  fkpj       !:
     2548       REAL(wp) ::  fkpjp      !:
     2549       REAL(wp) ::  logratio   !:
     2550       REAL(wp) ::  logzuc1    !:
     2551       REAL(wp) ::  zuc1       !:
     2552
     2553
     2554       ib = nxl
     2555       ie = nxr
     2556       jb = nys
     2557       je = nyn
     2558       IF ( nest_bound_l )  THEN
     2559          ib = nxl - 1
     2560!
     2561!--       For u, nxl is a ghost node, but not for the other variables
     2562          IF ( var == 'u' )  THEN
     2563             ib = nxl
     2564          ENDIF
     2565       ENDIF
     2566       IF ( nest_bound_s )  THEN
     2567          jb = nys - 1
     2568!
     2569!--       For v, nys is a ghost node, but not for the other variables
     2570          IF ( var == 'v' )  THEN
     2571             jb = nys
     2572          ENDIF
     2573       ENDIF
     2574       IF ( nest_bound_r )  THEN
     2575          ie = nxr + 1
     2576       ENDIF
     2577       IF ( nest_bound_n )  THEN
     2578          je = nyn + 1
     2579       ENDIF
     2580
     2581!
     2582!--    Trilinear interpolation.
     2583       DO  i = ib, ie
     2584          DO  j = jb, je
     2585             DO  k = kb(j,i), nzt + 1
     2586                l = ic(i)
     2587                m = jc(j)
     2588                n = kc(k)
     2589                fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
     2590                fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
     2591                fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
     2592                fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
     2593                fk       = r1y(j) * fkj  + r2y(j) * fkjp
     2594                fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
     2595                f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
     2596             ENDDO
     2597          ENDDO
     2598       ENDDO
     2599
     2600!
     2601!--    Correct the interpolated values of u and v in near-wall nodes, i.e. in
     2602!--    the nodes below the coarse-grid nodes with k=1. The corrction is only
     2603!--    made over horizontal wall surfaces in this phase. For the nest boundary
     2604!--    conditions, a corresponding correction is made for all vertical walls,
     2605!--    too.
     2606       IF ( var == 'u' .OR. var == 'v' )  THEN
     2607          DO  i = ib, nxr
     2608             DO  j = jb, nyn
     2609                kbc = 1
     2610!
     2611!--             kbc is the first coarse-grid point above the surface
     2612                DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )
     2613                   kbc = kbc + 1
     2614                ENDDO
     2615                zuc1 = cg%zu(kbc)
     2616                k1   = kb(j,i) + 1
     2617                DO  WHILE ( zu(k1) < zuc1 )
     2618                   k1 = k1 + 1
     2619                ENDDO
     2620                logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
     2621
     2622                k = kb(j,i) + 1
     2623                DO  WHILE ( zu(k) < zuc1 )
     2624                   logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) / logzuc1
     2625                   f(k,j,i) = logratio * f(k1,j,i)
     2626                   k  = k + 1
     2627                ENDDO
     2628                f(kb(j,i),j,i) = 0.0_wp
     2629             ENDDO
     2630          ENDDO
     2631
     2632       ELSEIF ( var == 'w' )  THEN
     2633
     2634          DO  i = ib, nxr
     2635              DO  j = jb, nyn
     2636                f(kb(j,i),j,i) = 0.0_wp
     2637             ENDDO
     2638          ENDDO
     2639
     2640       ENDIF
     2641
     2642    END SUBROUTINE pmci_interp_tril_all
     2643
     2644#endif
     2645 END SUBROUTINE pmci_client_initialize
     2646
     2647
     2648
     2649 SUBROUTINE pmci_ensure_nest_mass_conservation
     2650
     2651#if defined( __parallel )
     2652!
     2653!-- Adjust the volume-flow rate through the top boundary so that the net volume
     2654!-- flow through all boundaries of the current nest domain becomes zero.
     2655    IMPLICIT NONE
     2656
     2657    INTEGER(iwp) ::  i                          !:
     2658    INTEGER(iwp) ::  ierr                       !:
     2659    INTEGER(iwp) ::  j                          !:
     2660    INTEGER(iwp) ::  k                          !:
     2661
     2662    REAL(wp) ::  dxdy                            !:
     2663    REAL(wp) ::  innor                           !:
     2664    REAL(wp) ::  w_lt                            !:
     2665    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !:
     2666
     2667!
     2668!-- Sum up the volume flow through the left/right boundaries
     2669    volume_flow(1)   = 0.0_wp
     2670    volume_flow_l(1) = 0.0_wp
     2671
     2672    IF ( nest_bound_l )  THEN
     2673       i = 0
     2674       innor = dy
     2675       DO   j = nys, nyn
     2676          DO   k = nzb_u_inner(j,i)+1, nzt
     2677             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
     2678          ENDDO
     2679       ENDDO
     2680    ENDIF
     2681
     2682    IF ( nest_bound_r )  THEN
     2683       i = nx + 1
     2684       innor = -dy
     2685       DO   j = nys, nyn
     2686          DO   k = nzb_u_inner(j,i)+1, nzt
     2687             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)
     2688          ENDDO
     2689       ENDDO
     2690    ENDIF
     2691
     2692#if defined( __parallel )
     2693    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     2694    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
     2695                        MPI_SUM, comm2d, ierr )
     2696#else
     2697    volume_flow(1) = volume_flow_l(1)
     2698#endif
     2699
     2700!
     2701!-- Sum up the volume flow through the south/north boundaries
     2702    volume_flow(2)   = 0.0_wp
     2703    volume_flow_l(2) = 0.0_wp
     2704
     2705    IF ( nest_bound_s )  THEN
     2706       j = 0
     2707       innor = dx
     2708       DO   i = nxl, nxr
     2709          DO   k = nzb_v_inner(j,i)+1, nzt
     2710             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
     2711          ENDDO
     2712       ENDDO
     2713    ENDIF
     2714
     2715    IF ( nest_bound_n )  THEN
     2716       j = ny + 1
     2717       innor = -dx
     2718       DO   i = nxl, nxr
     2719          DO   k = nzb_v_inner(j,i)+1, nzt
     2720             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
     2721          ENDDO
     2722       ENDDO
     2723    ENDIF
     2724
     2725#if defined( __parallel )
     2726    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     2727    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL,         &
     2728                        MPI_SUM, comm2d, ierr )
     2729#else
     2730    volume_flow(2) = volume_flow_l(2)
     2731#endif
     2732
     2733!
     2734!-- Sum up the volume flow through the top boundary
     2735    volume_flow(3)   = 0.0_wp
     2736    volume_flow_l(3) = 0.0_wp
     2737    dxdy = dx * dy
     2738    k = nzt
     2739    DO   i = nxl, nxr
     2740       DO   j = nys, nyn
     2741          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
     2742       ENDDO
     2743    ENDDO
     2744
     2745#if defined( __parallel )
     2746    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     2747    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL,         &
     2748                        MPI_SUM, comm2d, ierr )
     2749#else
     2750    volume_flow(3) = volume_flow_l(3)
     2751#endif
     2752
     2753!
     2754!-- Correct the top-boundary value of w
     2755    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
     2756    DO   i = nxl, nxr
     2757       DO   j = nys, nyn
     2758          DO  k = nzt, nzt + 1
     2759             w(k,j,i) = w(k,j,i) + w_lt
     2760          ENDDO
     2761       ENDDO
     2762    ENDDO
     2763
     2764#endif
     2765 END SUBROUTINE pmci_ensure_nest_mass_conservation
     2766
     2767
     2768!-- TO_DO: the timestep sycnchronization could be done easier using
     2769!--        an MPI_ALLREDUCE with MIN over MPI_COMM_WORLD
    21752770 SUBROUTINE pmci_server_synchronize
    21762771
    21772772#if defined( __parallel )
    21782773!
    2179 !-- Unify the time steps for each model and synchronize.
    2180 !-- This is based on the assumption that the native time step
    2181 !-- (original dt_3d) of any server is always larger than the smallest
    2182 !-- native time step of it s clients.
     2774!-- Unify the time steps for each model and synchronize. This is based on the
     2775!-- assumption that the native time step (original dt_3d) of any server is
     2776!-- always larger than the smallest native time step of it s clients.
    21832777    IMPLICIT NONE
    2184     INTEGER(iwp)           ::  client_id   !:
     2778
     2779    INTEGER(iwp) ::  client_id   !:
     2780    INTEGER(iwp) ::  ierr        !:
     2781    INTEGER(iwp) ::  m           !:
     2782
    21852783    REAL(wp), DIMENSION(1) ::  dtc         !:
    21862784    REAL(wp), DIMENSION(1) ::  dtl         !:
    2187     INTEGER(iwp)           ::  ierr        !:
    2188     INTEGER(iwp)           ::  m           !:   
    2189 
    2190 !
    2191 !-- First find the smallest native time step of all the clients of the current server.
     2785
     2786!
     2787!-- First find the smallest native time step of all the clients of the current
     2788!-- server.
    21922789    dtl(1) = 999999.9_wp
    2193     DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
     2790    DO  m = 1, SIZE( PMC_Server_for_Client )-1
    21942791       client_id = PMC_Server_for_Client(m)
    21952792       IF ( myid == 0 )  THEN
     
    22012798
    22022799!
    2203 !-- Broadcast the unified time step to all server processes.
     2800!-- Broadcast the unified time step to all server processes
    22042801    CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
    22052802
    22062803!
    2207 !-- Send the new time step to all the clients of the current server.
     2804!-- Send the new time step to all the clients of the current server
    22082805    DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
    22092806       client_id = PMC_Server_for_Client(m)
     
    22222819#if defined( __parallel )
    22232820!
    2224 !-- Unify the time steps for each model and synchronize.
    2225 !-- This is based on the assumption that the native time step
    2226 !-- (original dt_3d) of any server is always larger than the smallest
    2227 !-- native time step of it s clients.
     2821!-- Unify the time steps for each model and synchronize. This is based on the
     2822!-- assumption that the native time step (original dt_3d) of any server is
     2823!-- always larger than the smallest native time step of it s clients.
    22282824
    22292825    IMPLICIT NONE
     2826
     2827    INTEGER(iwp) ::  ierr   !:
     2828
    22302829    REAL(wp), DIMENSION(1) ::  dtl    !:
    22312830    REAL(wp), DIMENSION(1) ::  dts    !:
    2232     INTEGER(iwp)           ::  ierr   !:
    22332831   
    22342832
    22352833    dtl(1) = dt_3d
    2236     IF ( cpl_id > 1 )  THEN   !  Root id is never a client
     2834    IF ( cpl_id > 1 )  THEN
    22372835       IF ( myid==0 )  THEN
    22382836          CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
     
    22422840
    22432841!
    2244 !--    Broadcast the unified time step to all server processes.
     2842!--    Broadcast the unified time step to all server processes
    22452843       CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
    22462844    ENDIF
     
    22552853    IMPLICIT NONE
    22562854
    2257     INTEGER(iwp),INTENT(IN) ::  swaplevel  !: swaplevel (1 or 2) of PALM's timestep
     2855    INTEGER(iwp),INTENT(IN) ::  swaplevel  !: swaplevel (1 or 2) of PALM's
     2856                                           !: timestep
    22582857
    22592858    INTEGER(iwp)            ::  client_id  !:
     
    22782877
    22792878#if defined( __parallel )
    2280     INTEGER(iwp)            ::  client_id   !:
    2281     INTEGER(iwp)            ::  i           !:
    2282     INTEGER(iwp)            ::  j           !:
    2283     INTEGER(iwp)            ::  ierr        !:
    2284     INTEGER(iwp)            ::  m           !:
    2285     REAL(wp)                ::  waittime    !:
    2286     REAL(wp), DIMENSION(1)  ::  dtc         !:
    2287     REAL(wp), DIMENSION(1)  ::  dtl         !:
    2288 
    2289 !
    2290 !-- First find the smallest native time step of all the clients of the current server.
     2879    INTEGER(iwp) ::  client_id   !:
     2880    INTEGER(iwp) ::  i           !:
     2881    INTEGER(iwp) ::  j           !:
     2882    INTEGER(iwp) ::  ierr        !:
     2883    INTEGER(iwp) ::  m           !:
     2884
     2885    REAL(wp)               ::  waittime    !:
     2886    REAL(wp), DIMENSION(1) ::  dtc         !:
     2887    REAL(wp), DIMENSION(1) ::  dtl         !:
     2888
     2889!
     2890!-- First find the smallest native time step of all the clients of the current
     2891!-- server.
    22912892    dtl(1) = 999999.9_wp
    2292     DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
    2293        client_id = PMC_Server_for_Client(m)
     2893    DO  m = 1, SIZE( pmc_server_for_client )-1
     2894       client_id = pmc_server_for_client(m)
    22942895       IF ( myid==0 )  THEN
    22952896          CALL pmc_recv_from_client( client_id, dtc, SIZE( dtc ), 0, 101, ierr )
     
    23002901
    23012902!
    2302 !-- Broadcast the unified time step to all server processes.
     2903!-- Broadcast the unified time step to all server processes
    23032904    CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
    23042905
    2305     DO  m = 1, SIZE( PMC_Server_for_Client ) - 1
     2906    DO  m = 1, SIZE( PMC_Server_for_Client )-1
    23062907       client_id = PMC_Server_for_Client(m)
    2307        CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
    2308 
    2309 !
    2310 !--    Send the new time step to all the clients of the current server.
    2311        IF ( myid == 0 ) THEN
     2908       CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )
     2909
     2910!
     2911!--    Send the new time step to all the clients of the current server
     2912       IF ( myid == 0 )  THEN
    23122913          CALL pmc_send_to_client( client_id, dtl, SIZE( dtl ), 0, 102, ierr )
    23132914       ENDIF
    2314        CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
     2915       CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )
    23152916       
    23162917       IF ( direction == server_to_client )  THEN
    2317           CALL cpu_log( log_point_s(71), 'PMC Server Send', 'start' )
     2918          CALL cpu_log( log_point_s(71), 'pmc server send', 'start' )
    23182919          CALL pmc_s_fillbuffer( client_id, waittime=waittime )
    2319           CALL cpu_log( log_point_s(71), 'PMC Server Send', 'stop' )
    2320        ELSE   !  Communication from client to server.
    2321           CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'start' )
     2920          CALL cpu_log( log_point_s(71), 'pmc server send', 'stop' )
     2921       ELSE
     2922!
     2923!--       Communication from client to server
     2924          CALL cpu_log( log_point_s(72), 'pmc server recv', 'start' )
    23222925          client_id = pmc_server_for_client(m)
    23232926          CALL pmc_s_getdata_from_buffer( client_id )
    2324           CALL cpu_log( log_point_s(72), 'PMC Server Recv', 'stop' )       
    2325 
    2326 !
    2327 !--       The anterpolated data is now available in u etc.
     2927          CALL cpu_log( log_point_s(72), 'pmc server recv', 'stop' )
     2928
     2929!
     2930!--       The anterpolated data is now available in u etc
    23282931          IF ( topography /= 'flat' )  THEN
    23292932
    23302933!
    23312934!--          Inside buildings/topography reset velocities and TKE back to zero.
    2332 !--          Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
    2333 !--          maybe revise later.
     2935!--          TO_DO: at least temperature should be included here immediately
     2936!--          Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at
     2937!--          present, maybe revise later.
    23342938             DO   i = nxlg, nxrg
    23352939                DO   j = nysg, nyng
     
    23372941                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
    23382942                   w(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
    2339                    e(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
     2943                   e(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
    23402944                ENDDO
    23412945             ENDDO
     
    23562960
    23572961#if defined( __parallel )
    2358     INTEGER(iwp)             ::  ierr        !:
    2359     INTEGER(iwp)             ::  icl         !:
    2360     INTEGER(iwp)             ::  icr         !:
    2361     INTEGER(iwp)             ::  jcs         !:
    2362     INTEGER(iwp)             ::  jcn         !:
     2962    INTEGER(iwp) ::  ierr        !:
     2963    INTEGER(iwp) ::  icl         !:
     2964    INTEGER(iwp) ::  icr         !:
     2965    INTEGER(iwp) ::  jcs         !:
     2966    INTEGER(iwp) ::  jcn         !:
    23632967   
    2364     REAL(wp), DIMENSION(1)   ::  dtl         !:
    2365     REAL(wp), DIMENSION(1)   ::  dts         !:
    2366     REAL(wp)                 ::  waittime    !:
     2968    REAL(wp) ::  waittime    !:
     2969
     2970    REAL(wp), DIMENSION(1) ::  dtl         !:
     2971    REAL(wp), DIMENSION(1) ::  dts         !:
    23672972
    23682973
    23692974    dtl = dt_3d
    2370     IF ( cpl_id > 1 )  THEN   !  Root id is never a client     
    2371        CALL cpu_log( log_point_s(70), 'PMC model sync', 'start' )
     2975    IF ( cpl_id > 1 )  THEN
     2976       CALL cpu_log( log_point_s(70), 'pmc sync', 'start' )
    23722977       IF ( myid==0 )  THEN
    23732978          CALL pmc_send_to_server( dtl, SIZE( dtl ), 0, 101, ierr )
     
    23792984!--    Broadcast the unified time step to all server processes.
    23802985       CALL MPI_BCAST( dt_3d, 1, MPI_REAL, 0, comm2d, ierr )
    2381        CALL cpu_log( log_point_s(70), 'PMC model sync', 'stop' )
     2986       CALL cpu_log( log_point_s(70), 'pmc sync', 'stop' )
    23822987
    23832988!
     
    23892994
    23902995       IF ( direction == server_to_client )  THEN
    2391           CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'start' )
     2996
     2997          CALL cpu_log( log_point_s(73), 'pmc client recv', 'start' )
    23922998          CALL pmc_c_getbuffer( WaitTime = WaitTime )
    2393           CALL cpu_log( log_point_s(73), 'PMC Client Recv', 'stop' )
     2999          CALL cpu_log( log_point_s(73), 'pmc client recv', 'stop' )
    23943000
    23953001          CALL pmci_interpolation
    23963002
    2397        ELSE    !  IF ( direction == server_to_client )
    2398 
     3003       ELSE
     3004!
     3005!--       direction == server_to_client
    23993006          CALL pmci_anterpolation
    24003007
    2401           CALL cpu_log( log_point_s(74), 'PMC Client Send', 'start' )
     3008          CALL cpu_log( log_point_s(74), 'pmc client send', 'start' )
    24023009          CALL pmc_c_putbuffer( WaitTime = WaitTime )
    2403           CALL cpu_log( log_point_s(74), 'PMC Client Send', 'stop' )
     3010          CALL cpu_log( log_point_s(74), 'pmc client send', 'stop' )
     3011
    24043012       ENDIF
    24053013    ENDIF
     
    24073015 CONTAINS
    24083016
    2409 
    2410    SUBROUTINE pmci_interpolation
    2411 
    2412 !
    2413 !--   A wrapper routine for all interpolation and extrapolation actions.
    2414       IMPLICIT NONE
    2415 
    2416 !
    2417 !--   Add IF-condition here: IF not vertical nesting
    2418       IF ( nest_bound_l )  THEN   !  Left border pe
    2419          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
    2420                                    nzb_u_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'u' )
    2421          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
    2422                                    nzb_v_inner, logc_v_l, logc_ratio_v_l, nzt_topo_nestbc_l, 'l', 'v' )
    2423          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
    2424                                    nzb_w_inner, logc_w_l, logc_ratio_w_l, nzt_topo_nestbc_l, 'l', 'w' )
    2425          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2426                                    nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 'e' )
    2427          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2428                                    nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
    2429          IF ( humidity  .OR.  passive_scalar )  THEN
    2430             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
    2431                                       nzb_s_inner, logc_u_l, logc_ratio_u_l, nzt_topo_nestbc_l, 'l', 's' )
    2432          ENDIF
    2433          IF ( nesting_mode == 'one-way' )  THEN
    2434             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
    2435             CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
    2436             CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
    2437             CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
    2438             CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
    2439             IF ( humidity  .OR.  passive_scalar )  THEN
    2440                CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
    2441             ENDIF
    2442          ENDIF
    2443       ENDIF
    2444       IF ( nest_bound_r )  THEN   !  Right border pe
    2445          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
    2446                                    nzb_u_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'u' )
    2447          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
    2448                                    nzb_v_inner, logc_v_r, logc_ratio_v_r, nzt_topo_nestbc_r, 'r', 'v' )
    2449          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
    2450                                    nzb_w_inner, logc_w_r, logc_ratio_w_r, nzt_topo_nestbc_r, 'r', 'w' )
    2451          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2452                                    nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 'e' )
    2453          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2454                                    nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
    2455          IF ( humidity  .OR.  passive_scalar )  THEN
    2456             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
    2457                                       nzb_s_inner, logc_u_r, logc_ratio_u_r, nzt_topo_nestbc_r, 'r', 's' )
    2458          ENDIF
    2459          IF ( nesting_mode == 'one-way' )  THEN
    2460             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
     3017    SUBROUTINE pmci_interpolation
     3018
     3019!
     3020!--    A wrapper routine for all interpolation and extrapolation actions
     3021       IMPLICIT NONE
     3022
     3023!
     3024!--    Add IF-condition here: IF not vertical nesting
     3025!--    Left border pe:
     3026       IF ( nest_bound_l )  THEN
     3027          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
     3028                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_l,   &
     3029                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
     3030                                    'u' )
     3031          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
     3032                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_l,   &
     3033                                    logc_ratio_v_l, nzt_topo_nestbc_l, 'l',    &
     3034                                    'v' )
     3035          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
     3036                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_l,   &
     3037                                    logc_ratio_w_l, nzt_topo_nestbc_l, 'l',    &
     3038                                    'w' )
     3039          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
     3040                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_l,   &
     3041                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
     3042                                    'e' )
     3043          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
     3044                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_l,   &
     3045                                    logc_ratio_u_l, nzt_topo_nestbc_l, 'l',    &
     3046                                    's' )
     3047          IF ( humidity  .OR.  passive_scalar )  THEN
     3048             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
     3049                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
     3050                                       logc_u_l, logc_ratio_u_l,               &
     3051                                       nzt_topo_nestbc_l, 'l', 's' )
     3052          ENDIF
     3053
     3054          IF ( nesting_mode == 'one-way' )  THEN
     3055             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'l', 'u' )
     3056             CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'l', 'v' )
     3057             CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'l', 'w' )
     3058             CALL pmci_extrap_ifoutflow_lr( e, nzb_s_inner, 'l', 'e' )
     3059             CALL pmci_extrap_ifoutflow_lr( pt,nzb_s_inner, 'l', 's' )
     3060             IF ( humidity  .OR.  passive_scalar )  THEN
     3061                CALL pmci_extrap_ifoutflow_lr( q, nzb_s_inner, 'l', 's' )
     3062             ENDIF
     3063          ENDIF
     3064
     3065       ENDIF
     3066!
     3067!--    Right border pe
     3068       IF ( nest_bound_r )  THEN
     3069          CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
     3070                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_r,   &
     3071                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
     3072                                    'u' )
     3073          CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
     3074                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_r,   &
     3075                                    logc_ratio_v_r, nzt_topo_nestbc_r, 'r',    &
     3076                                    'v' )
     3077          CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
     3078                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_r,   &
     3079                                    logc_ratio_w_r, nzt_topo_nestbc_r, 'r',    &
     3080                                    'w' )
     3081          CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
     3082                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_r,   &
     3083                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
     3084                                    'e' )
     3085          CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
     3086                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_r,   &
     3087                                    logc_ratio_u_r, nzt_topo_nestbc_r, 'r',    &
     3088                                    's' )
     3089          IF ( humidity  .OR.  passive_scalar )  THEN
     3090             CALL pmci_interp_tril_lr( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
     3091                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
     3092                                       logc_u_r, logc_ratio_u_r,               &
     3093                                       nzt_topo_nestbc_r, 'r', 's' )
     3094          ENDIF
     3095
     3096          IF ( nesting_mode == 'one-way' )  THEN
     3097             CALL pmci_extrap_ifoutflow_lr( u, nzb_u_inner, 'r', 'u' )
    24613098             CALL pmci_extrap_ifoutflow_lr( v, nzb_v_inner, 'r', 'v' )
    24623099             CALL pmci_extrap_ifoutflow_lr( w, nzb_w_inner, 'r', 'w' )
     
    24673104             ENDIF
    24683105          ENDIF
    2469        ENDIF
    2470        IF ( nest_bound_s )  THEN   !  South border pe
    2471           CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
    2472                                     nzb_u_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'u' )
    2473           CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
    2474                                     nzb_v_inner, logc_v_s, logc_ratio_v_s, nzt_topo_nestbc_s, 's', 'v' )
    2475           CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
    2476                                     nzb_w_inner, logc_w_s, logc_ratio_w_s, nzt_topo_nestbc_s, 's', 'w' )
    2477           CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2478                                     nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 'e' )
    2479           CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2480                                     nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
     3106
     3107       ENDIF
     3108!
     3109!--    South border pe
     3110       IF ( nest_bound_s )  THEN
     3111          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
     3112                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_s,   &
     3113                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
     3114                                    'u' )
     3115          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
     3116                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_s,   &
     3117                                    logc_ratio_v_s, nzt_topo_nestbc_s, 's',    &
     3118                                    'v' )
     3119          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
     3120                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_s,   &
     3121                                    logc_ratio_w_s, nzt_topo_nestbc_s, 's',    &
     3122                                    'w' )
     3123          CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
     3124                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_s,   &
     3125                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
     3126                                    'e' )
     3127          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
     3128                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_s,   &
     3129                                    logc_ratio_u_s, nzt_topo_nestbc_s, 's',    &
     3130                                    's' )
    24813131          IF ( humidity  .OR.  passive_scalar )  THEN
    2482              CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,         &
    2483                                        nzb_s_inner, logc_u_s, logc_ratio_u_s, nzt_topo_nestbc_s, 's', 's' )
     3132             CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
     3133                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
     3134                                       logc_u_s, logc_ratio_u_s,               &
     3135                                       nzt_topo_nestbc_s, 's', 's' )
    24843136          ENDIF
     3137
    24853138          IF ( nesting_mode == 'one-way' )  THEN
    24863139             CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 's', 'u' )
     
    24933146             ENDIF
    24943147          ENDIF
    2495        ENDIF
    2496        IF ( nest_bound_n )  THEN   !  North border pe
    2497           CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo,          &
    2498                                     nzb_u_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'u' )
    2499           CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo,          &
    2500                                     nzb_v_inner, logc_v_n, logc_ratio_v_n, nzt_topo_nestbc_n, 'n', 'v' )
    2501           CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw,          &
    2502                                     nzb_w_inner, logc_w_n, logc_ratio_w_n, nzt_topo_nestbc_n, 'n', 'w' )
    2503           CALL pmci_interp_tril_sn( e,  ec,  ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                  &
    2504                                     nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 'e' )
    2505           CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo,          &
    2506                                     nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
     3148
     3149       ENDIF
     3150!
     3151!--    North border pe
     3152       IF ( nest_bound_n )  THEN
     3153          CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,  &
     3154                                    r2yo, r1zo, r2zo, nzb_u_inner, logc_u_n,   &
     3155                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
     3156                                    'u' )
     3157          CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,  &
     3158                                    r2yv, r1zo, r2zo, nzb_v_inner, logc_v_n,   &
     3159                                    logc_ratio_v_n, nzt_topo_nestbc_n, 'n',    &
     3160                                    'v' )
     3161          CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,  &
     3162                                    r2yo, r1zw, r2zw, nzb_w_inner, logc_w_n,   &
     3163                                    logc_ratio_w_n, nzt_topo_nestbc_n, 'n',    &
     3164                                    'w' )
     3165          CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,  &
     3166                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_n,   &
     3167                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
     3168                                    'e' )
     3169          CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,  &
     3170                                    r2yo, r1zo, r2zo, nzb_s_inner, logc_u_n,   &
     3171                                    logc_ratio_u_n, nzt_topo_nestbc_n, 'n',    &
     3172                                    's' )
    25073173          IF ( humidity  .OR.  passive_scalar )  THEN
    2508              CALL pmci_interp_tril_sn( q, qc, ico,jco,kco,r1xo,r2xo,r1yo,r2yo,r1zo,r2zo,                 &
    2509                                        nzb_s_inner, logc_u_n, logc_ratio_u_n, nzt_topo_nestbc_n, 'n', 's' )
     3174             CALL pmci_interp_tril_sn( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, &
     3175                                       r2yo, r1zo, r2zo, nzb_s_inner,          &
     3176                                       logc_u_n, logc_ratio_u_n,               &
     3177                                       nzt_topo_nestbc_n, 'n', 's' )
    25103178          ENDIF
     3179
    25113180          IF ( nesting_mode == 'one-way' )  THEN
    25123181             CALL pmci_extrap_ifoutflow_sn( u, nzb_u_inner, 'n', 'u' )
     
    25183187                CALL pmci_extrap_ifoutflow_sn( q, nzb_s_inner, 'n', 's' )
    25193188             ENDIF
     3189
    25203190          ENDIF
    2521       ENDIF
    2522 
    2523 !
    2524 !--   All PEs are top-border PEs
    2525        CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, 'u' )
    2526        CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, 'v' )
    2527        CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, 'w' )
    2528        CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 'e' )
    2529        CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
    2530        IF ( humidity  .OR.  passive_scalar )  THEN
    2531           CALL pmci_interp_tril_t( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, 's' )
    2532        ENDIF
     3191
     3192       ENDIF
     3193
     3194!
     3195!--    All PEs are top-border PEs
     3196       CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
     3197                                r2yo, r1zo, r2zo, 'u' )
     3198       CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,      &
     3199                                r2yv, r1zo, r2zo, 'v' )
     3200       CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
     3201                                r2yo, r1zw, r2zw, 'w' )
     3202       CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,      &
     3203                                r2yo, r1zo, r2zo, 'e' )
     3204       CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,      &
     3205                                r2yo, r1zo, r2zo, 's' )
     3206       IF ( humidity .OR. passive_scalar )  THEN
     3207          CALL pmci_interp_tril_t( q, qc, ico, jco, kco, r1xo, r2xo, r1yo,     &
     3208                                   r2yo, r1zo, r2zo, 's' )
     3209       ENDIF
     3210
    25333211       IF ( nesting_mode == 'one-way' )  THEN
    25343212          CALL pmci_extrap_ifoutflow_t( u,  'u' )
     
    25413219          ENDIF
    25423220      ENDIF
     3221
    25433222   END SUBROUTINE pmci_interpolation
    25443223
     
    25513230      IMPLICIT NONE
    25523231
    2553       CALL pmci_anterp_tophat( u,  uc,  kceu, iflu, ifuu, jflo, jfuo, kflo, kfuo, nzb_u_inner, 'u' )
    2554       CALL pmci_anterp_tophat( v,  vc,  kceu, iflo, ifuo, jflv, jfuv, kflo, kfuo, nzb_v_inner, 'v' )
    2555       CALL pmci_anterp_tophat( w,  wc,  kcew, iflo, ifuo, jflo, jfuo, kflw, kfuw, nzb_w_inner, 'w' )
    2556       CALL pmci_anterp_tophat( pt, ptc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
     3232      CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
     3233                               kfuo, nzb_u_inner, 'u' )
     3234      CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
     3235                               kfuo, nzb_v_inner, 'v' )
     3236      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
     3237                               kfuw, nzb_w_inner, 'w' )
     3238      CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo,    &
     3239                               kfuo, nzb_s_inner, 's' )
    25573240      IF ( humidity  .OR.  passive_scalar )  THEN
    2558          CALL pmci_anterp_tophat( q, qc, kceu, iflo, ifuo, jflo, jfuo, kflo, kfuo, nzb_s_inner, 's' )
     3241         CALL pmci_anterp_tophat( q, qc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
     3242                                  kfuo, nzb_s_inner, 's' )
    25593243      ENDIF
    25603244
     
    25633247
    25643248
    2565    SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, logc, logc_ratio, &
    2566                                   nzt_topo_nestbc, edge, var )
    2567 
     3249   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
     3250                                   r2z, kb, logc, logc_ratio, nzt_topo_nestbc, &
     3251                                   edge, var )
    25683252!
    25693253!--   Interpolation of ghost-node values used as the client-domain boundary
    2570 !--   conditions. This subroutine handles the left and right boundaries.
    2571 !--   This subroutine is based on trilinear interpolation.
    2572 !--   Constant dz is still assumed.
    2573 !
    2574 !--                                   Antti Hellsten 22.2.2015.
    2575 !
    2576 !--   Rewritten so that all the coefficients and client-array indices are
    2577 !--   precomputed in the initialization phase by pmci_init_interp_tril.
    2578 !
    2579 !                                     Antti Hellsten 3.3.2015.
    2580 !
    2581 !--   Constant dz no more assumed.
    2582 !                                     Antti Hellsten 23.3.2015.
    2583 !
    2584 !--   Adapted for non-flat topography. However, the near-wall velocities
    2585 !--   are log-corrected only over horizontal surfaces, not yet near vertical
    2586 !--   walls.
    2587 !--                                   Antti Hellsten 26.3.2015.
    2588 !
    2589 !--   Indexing in the principal direction (i) is changed. Now, the nest-boundary
    2590 !--   values are interpolated only into the first ghost-node layers on each later
    2591 !--   boundary. These values are then simply copied to the second ghost-node layer.
    2592 !
    2593 !--                                   Antti Hellsten 6.10.2015.
     3254!--   conditions. This subroutine handles the left and right boundaries. It is
     3255!--   based on trilinear interpolation. Constant dz is still assumed.
     3256!--   TO_DO:  constant dz is an important restriction and should be checked
     3257!--           somewhere in order to let users not run into this trap
    25943258      IMPLICIT NONE
     3259
     3260!--   TO_DO: wrap long lines in this and the remaining interp_tril routines
    25953261      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT)          ::  f            !:
    25963262      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)                 ::  fc           !:
     
    26153281      INTEGER(iwp) ::  i       !:
    26163282      INTEGER(iwp) ::  ib      !:
     3283      INTEGER(iwp) ::  ibgp    !:
    26173284      INTEGER(iwp) ::  iw      !:
    26183285      INTEGER(iwp) ::  j       !:
     
    27593426
    27603427!
    2761 !--   Store the boundary values also into the second ghost node layer.
    2762       f(0:nzt+1,nys:nyn+1,ib) = f(0:nzt+1,nys:nyn+1,i)
     3428!--   Store the boundary values also into the other redundant ghost node layers
     3429      IF ( edge == 'l' )  THEN
     3430         DO ibgp = -nbgp, ib
     3431            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
     3432         ENDDO
     3433      ELSE IF ( edge == 'r' )  THEN
     3434         DO ibgp = ib, nx + nbgp
     3435            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
     3436         ENDDO
     3437      ENDIF
    27633438
    27643439   END SUBROUTINE pmci_interp_tril_lr
     
    28233498      INTEGER(iwp) ::  j       !:
    28243499      INTEGER(iwp) ::  jb      !:
     3500      INTEGER(iwp) ::  jbgp    !:
    28253501      INTEGER(iwp) ::  k       !:
    28263502      INTEGER(iwp) ::  kcorr   !:
     
    29603636
    29613637!
    2962 !--   Store the boundary values also into the second ghost node layer.
    2963       f(0:nzt+1,jb,nxl:nxr+1) = f(0:nzt+1,j,nxl:nxr+1)
     3638!--   Store the boundary values also into the other redundant ghost node layers
     3639      IF ( edge == 's' )  THEN
     3640         DO jbgp = -nbgp, jb
     3641            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
     3642         ENDDO
     3643      ELSE IF ( edge == 'n' )  THEN
     3644         DO jbgp = jb, ny + nbgp
     3645            f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
     3646         ENDDO
     3647      ENDIF
    29643648
    29653649   END SUBROUTINE pmci_interp_tril_sn
     
    29743658!--   This subroutine is based on trilinear interpolation.
    29753659!--   Constant dz is still assumed.
    2976 !
    2977 !--                                      Antti Hellsten 23.2.2015.
    2978 !
    2979 !
    2980 !--   Rewritten so that all the coefficients and client-array indices are
    2981 !--   precomputed in the initialization phase by pmci_init_interp_tril.
    2982 !
    2983 !--                                      Antti Hellsten 3.3.2015.
    2984 !
    2985 !--   Constant dz no more assumed.
    2986 !--                                      Antti Hellsten 23.3.2015.
    2987 !
    2988 !--   Indexing in the principal direction (k) is changed. Now, the nest-boundary
    2989 !--   values are interpolated only into the first ghost-node layer. Actually there is
    2990 !--   only one ghost-node layer in the k-direction. 
    2991 !
    2992 !--                                      Antti Hellsten 6.10.2015.
    2993 !
    29943660      IMPLICIT NONE
    29953661      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
     
    30693735
    30703736
    3071    SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
    3072 
    3073 !
    3074 !--   After the interpolation of ghost-node values for the client-domain boundary
    3075 !--   conditions, this subroutine checks if there is a local outflow through the
    3076 !--   boundary. In that case this subroutine overwrites the interpolated values
    3077 !--   by values extrapolated from the domain. This subroutine handles the left and
    3078 !--   right boundaries.
    3079 !--   However, this operation is only needed in case of one-way coupling.
    3080 !
    3081 !--                                       Antti Hellsten 9.3.2015.
    3082 !
    3083 !--   Indexing in the principal direction (i) is changed. Now, the nest-boundary
    3084 !--   values are interpolated only into the first ghost-node layers on each later
    3085 !--   boundary. These values are then simply copied to the second ghost-node layer.
    3086 !
    3087 !--                                       Antti Hellsten 6.10.2015.
    3088 !
    3089       IMPLICIT NONE
    3090       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
    3091       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
    3092 
    3093       CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
    3094       CHARACTER(LEN=1),INTENT(IN) ::  var    !:
    3095 
    3096       INTEGER(iwp) ::  i     !:
    3097       INTEGER(iwp) ::  ib    !:
    3098       INTEGER(iwp) ::  ied   !:
    3099       INTEGER(iwp) ::  j     !:
    3100       INTEGER(iwp) ::  k     !:
     3737    SUBROUTINE pmci_extrap_ifoutflow_lr( f, kb, edge, var )
     3738!
     3739!--    After the interpolation of ghost-node values for the client-domain
     3740!--    boundary conditions, this subroutine checks if there is a local outflow
     3741!--    through the boundary. In that case this subroutine overwrites the
     3742!--    interpolated values by values extrapolated from the domain. This
     3743!--    subroutine handles the left and right boundaries. However, this operation
     3744!--    is only needed in case of one-way coupling.
     3745       IMPLICIT NONE
     3746
     3747       CHARACTER(LEN=1),INTENT(IN) ::  edge   !:
     3748       CHARACTER(LEN=1),INTENT(IN) ::  var    !:
     3749
     3750       INTEGER(iwp) ::  i     !:
     3751       INTEGER(iwp) ::  ib    !:
     3752       INTEGER(iwp) ::  ibgp  !:
     3753       INTEGER(iwp) ::  ied   !:
     3754       INTEGER(iwp) ::  j     !:
     3755       INTEGER(iwp) ::  k     !:
    31013756     
    3102       REAL(wp) ::  outnor    !:
    3103       REAL(wp) ::  vdotnor   !:
    3104 
    3105 !
    3106 !--   Check which edge is to be handled: left or right.
    3107       IF ( edge == 'l' )  THEN
    3108          IF ( var == 'u' )  THEN
    3109             i   = nxl
    3110             ib  = nxl - 1
    3111             ied = nxl + 1
    3112          ELSE
    3113             i   = nxl - 1
    3114             ib  = nxl - 2
    3115             ied = nxl
    3116          ENDIF
    3117          outnor = -1.0_wp
    3118       ELSEIF ( edge == 'r' )  THEN
    3119          i      = nxr + 1
    3120          ib     = nxr + 2
    3121          ied    = nxr
    3122          outnor = 1.0_wp
    3123       ENDIF
    3124 
    3125       DO  j = nys, nyn + 1
    3126          DO  k = kb(j,i), nzt +1
    3127             vdotnor = outnor * u(k,j,ied)
    3128             IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
    3129                f(k,j,i) = f(k,j,ied)
    3130             ENDIF
    3131          ENDDO
    3132          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
    3133             f(kb(j,i),j,i) = 0.0_wp
    3134          ENDIF
    3135       ENDDO
    3136 
    3137 !
    3138 !--   Store the updated boundary values also into the second ghost node layer.
    3139       f(0:nzt,nys:nyn+1,ib) = f(0:nzt,nys:nyn+1,i)
    3140 
    3141    END SUBROUTINE pmci_extrap_ifoutflow_lr
    3142 
    3143 
    3144 
    3145    SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
    3146 !
    3147 !--   After  the interpolation of ghost-node values for the client-domain boundary
    3148 !--   conditions, this subroutine checks if there is a local outflow through the
    3149 !--   boundary. In that case this subroutine overwrites the interpolated values
    3150 !--   by values extrapolated from the domain. This subroutine handles the south and
    3151 !--   north boundaries.
    3152 !
    3153 !--                                       Antti Hellsten 9.3.2015.
    3154 !
    3155 !--   Indexing in the principal direction (j) is changed. Now, the nest-boundary
    3156 !--   values are interpolated only into the first ghost-node layers on each later
    3157 !--   boundary. These values are then simply copied to the second ghost-node layer.
    3158 !
    3159 !--                                       Antti Hellsten 6.10.2015.
    3160 
    3161       IMPLICIT NONE
    3162       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f    !:
    3163       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb   !:
    3164       CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
    3165       CHARACTER(LEN=1), INTENT(IN) ::  var    !:
     3757       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
     3758
     3759       REAL(wp) ::  outnor    !:
     3760       REAL(wp) ::  vdotnor   !:
     3761
     3762       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
     3763
     3764!
     3765!--    Check which edge is to be handled: left or right
     3766       IF ( edge == 'l' )  THEN
     3767          IF ( var == 'u' )  THEN
     3768             i   = nxl
     3769             ib  = nxl - 1
     3770             ied = nxl + 1
     3771          ELSE
     3772             i   = nxl - 1
     3773             ib  = nxl - 2
     3774             ied = nxl
     3775          ENDIF
     3776          outnor = -1.0_wp
     3777       ELSEIF ( edge == 'r' )  THEN
     3778          i      = nxr + 1
     3779          ib     = nxr + 2
     3780          ied    = nxr
     3781          outnor = 1.0_wp
     3782       ENDIF
     3783
     3784       DO  j = nys, nyn + 1
     3785          DO  k = kb(j,i), nzt +1
     3786             vdotnor = outnor * u(k,j,ied)
     3787!
     3788!--          Local outflow
     3789             IF ( vdotnor > 0.0_wp )  THEN
     3790                f(k,j,i) = f(k,j,ied)
     3791             ENDIF
     3792          ENDDO
     3793          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
     3794             f(kb(j,i),j,i) = 0.0_wp
     3795          ENDIF
     3796       ENDDO
     3797
     3798!
     3799!--    Store the boundary values also into the redundant ghost node layers.
     3800       IF ( edge == 'l' )  THEN
     3801          DO ibgp = -nbgp, ib
     3802             f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
     3803          ENDDO
     3804       ELSEIF ( edge == 'r' )  THEN
     3805          DO ibgp = ib, nx + nbgp
     3806             f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
     3807          ENDDO
     3808       ENDIF
     3809
     3810    END SUBROUTINE pmci_extrap_ifoutflow_lr
     3811
     3812
     3813
     3814    SUBROUTINE pmci_extrap_ifoutflow_sn( f, kb, edge, var )
     3815!
     3816!--    After  the interpolation of ghost-node values for the client-domain
     3817!--    boundary conditions, this subroutine checks if there is a local outflow
     3818!--    through the boundary. In that case this subroutine overwrites the
     3819!--    interpolated values by values extrapolated from the domain. This
     3820!--    subroutine handles the south and north boundaries.
     3821       IMPLICIT NONE
     3822
     3823       CHARACTER(LEN=1), INTENT(IN) ::  edge   !:
     3824       CHARACTER(LEN=1), INTENT(IN) ::  var    !:
    31663825     
    3167       INTEGER(iwp) ::  i         !:
    3168       INTEGER(iwp) ::  j         !:
    3169       INTEGER(iwp) ::  jb        !:
    3170       INTEGER(iwp) ::  jed       !:
    3171       INTEGER(iwp) ::  k         !:
    3172       REAL(wp)     ::  outnor    !:
    3173       REAL(wp)     ::  vdotnor   !:
    3174 
    3175 !
    3176 !--   Check which edge is to be handled: left or right.
    3177       IF ( edge == 's' )  THEN
    3178          IF ( var == 'v' )  THEN
    3179             j   = nys
    3180             jb  = nys - 1
    3181             jed = nys + 1
    3182          ELSE
    3183             j   = nys - 1
    3184             jb  = nys - 2
    3185             jed = nys
    3186          ENDIF
    3187          outnor = -1.0_wp
    3188       ELSEIF ( edge == 'n' )  THEN
    3189          j      = nyn + 1
    3190          jb     = nyn + 2
    3191          jed    = nyn
    3192          outnor = 1.0_wp
    3193       ENDIF
    3194 
    3195       DO  i = nxl, nxr + 1
    3196          DO  k = kb(j,i), nzt + 1
    3197             vdotnor = outnor * v(k,jed,i)
    3198             IF ( vdotnor > 0.0_wp )  THEN   !  Local outflow.
    3199                f(k,j,i) = f(k,jed,i)                 
    3200             ENDIF
    3201          ENDDO
    3202          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
    3203             f(kb(j,i),j,i) = 0.0_wp
    3204          ENDIF
    3205       ENDDO
    3206 
    3207 !
    3208 !--   Store the updated boundary values also into the second ghost node layer.
    3209       f(0:nzt,jb,nxl:nxr+1) = f(0:nzt,j,nxl:nxr+1)
    3210 
    3211    END SUBROUTINE pmci_extrap_ifoutflow_sn
     3826       INTEGER(iwp) ::  i         !:
     3827       INTEGER(iwp) ::  j         !:
     3828       INTEGER(iwp) ::  jb        !:
     3829       INTEGER(iwp) ::  jbgp      !:
     3830       INTEGER(iwp) ::  jed       !:
     3831       INTEGER(iwp) ::  k         !:
     3832
     3833       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb   !:
     3834
     3835       REAL(wp)     ::  outnor    !:
     3836       REAL(wp)     ::  vdotnor   !:
     3837
     3838       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !:
     3839
     3840!
     3841!--    Check which edge is to be handled: left or right
     3842       IF ( edge == 's' )  THEN
     3843          IF ( var == 'v' )  THEN
     3844             j   = nys
     3845             jb  = nys - 1
     3846             jed = nys + 1
     3847          ELSE
     3848             j   = nys - 1
     3849             jb  = nys - 2
     3850             jed = nys
     3851          ENDIF
     3852          outnor = -1.0_wp
     3853       ELSEIF ( edge == 'n' )  THEN
     3854          j      = nyn + 1
     3855          jb     = nyn + 2
     3856          jed    = nyn
     3857          outnor = 1.0_wp
     3858       ENDIF
     3859
     3860       DO  i = nxl, nxr + 1
     3861          DO  k = kb(j,i), nzt + 1
     3862             vdotnor = outnor * v(k,jed,i)
     3863!
     3864!--          Local outflow
     3865             IF ( vdotnor > 0.0_wp )  THEN
     3866                f(k,j,i) = f(k,jed,i)
     3867             ENDIF
     3868          ENDDO
     3869          IF ( (var == 'u' )  .OR.  (var == 'v' )  .OR.  (var == 'w') )  THEN
     3870             f(kb(j,i),j,i) = 0.0_wp
     3871          ENDIF
     3872       ENDDO
     3873
     3874!
     3875!--    Store the boundary values also into the redundant ghost node layers.
     3876       IF ( edge == 's' )  THEN
     3877          DO jbgp = -nbgp, jb
     3878             f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
     3879          ENDDO
     3880       ELSEIF ( edge == 'n' )  THEN
     3881          DO jbgp = jb, ny + nbgp
     3882             f(0:nzt+1,jbgp,nxlg:nxrg) = f(0:nzt+1,j,nxlg:nxrg)
     3883          ENDDO
     3884       ENDIF
     3885
     3886    END SUBROUTINE pmci_extrap_ifoutflow_sn
    32123887
    32133888 
    32143889
    3215    SUBROUTINE pmci_extrap_ifoutflow_t( f, var )
    3216 
    3217 !
    3218 !--   Interpolation of ghost-node values used as the client-domain boundary
    3219 !--   conditions. This subroutine handles the top boundary.
    3220 !--   This subroutine is based on trilinear interpolation.
    3221 !     
    3222 !--                                       Antti Hellsten 23.2.2015.
    3223 !     
    3224 !     
    3225 !--   Rewritten so that all the coefficients and client-array indices are
    3226 !--   precomputed in the initialization phase by init_interp_tril.
    3227 !     
    3228 !--                                       Antti Hellsten 3.3.2015.
    3229 !     
    3230 !--   Indexing in the principal direction (k) is changed. Now, the nest-boundary
    3231 !--   values are extrapolated only into the first ghost-node layer. Actually there is
    3232 !--   only one ghost-node layer in the k-direction. 
    3233 !     
    3234 !--                                       Antti Hellsten 6.10.2015.
    3235       IMPLICIT NONE
    3236       REAL(wp), DIMENSION(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp), INTENT(INOUT) ::  f   !:
    3237       CHARACTER(LEN=1), INTENT(IN) ::  var   !:
     3890    SUBROUTINE pmci_extrap_ifoutflow_t( f, var )
     3891!
     3892!--    Interpolation of ghost-node values used as the client-domain boundary
     3893!--    conditions. This subroutine handles the top boundary. It is based on
     3894!--    trilinear interpolation.
     3895       IMPLICIT NONE
     3896
     3897       CHARACTER(LEN=1), INTENT(IN) ::  var   !:
    32383898     
    3239       INTEGER(iwp) ::  i     !:
    3240       INTEGER(iwp) ::  j     !:
    3241       INTEGER(iwp) ::  k     !:
    3242       INTEGER(iwp) ::  ked   !:
    3243       REAL(wp) ::  vdotnor   !:
     3899       INTEGER(iwp) ::  i     !:
     3900       INTEGER(iwp) ::  j     !:
     3901       INTEGER(iwp) ::  k     !:
     3902       INTEGER(iwp) ::  ked   !:
     3903
     3904       REAL(wp) ::  vdotnor   !:
     3905
     3906       REAL(wp), DIMENSION(nzb:nzt+1,nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp),     &
     3907                 INTENT(INOUT) ::  f   !:
    32443908     
    32453909
    3246       IF ( var == 'w' )  THEN
    3247          k    = nzt
    3248          ked  = nzt - 1
    3249       ELSE
    3250          k    = nzt + 1 
    3251          ked  = nzt
    3252       ENDIF
    3253 
    3254       DO  i = nxl, nxr
    3255          DO  j = nys, nyn
    3256             vdotnor = w(ked,j,i)
    3257             IF ( vdotnor > 0.0_wp )  THEN   !:  Local outflow.             
    3258                f(k,j,i) = f(ked,j,i)
    3259             ENDIF
    3260          ENDDO
    3261       ENDDO
    3262 
    3263 !
    3264 !--   Just fill up the second ghost-node layer for w.
    3265       IF ( var == 'w' )  THEN
    3266          f(nzt+1,:,:) = f(nzt,:,:)
    3267       ENDIF
    3268 
    3269    END SUBROUTINE pmci_extrap_ifoutflow_t
    3270 
    3271 
    3272 
    3273    SUBROUTINE pmci_anterp_tophat( f, fc, kce, ifl, ifu, jfl, jfu, kfl, kfu, kb, var )
    3274 !
    3275 !--   Anterpolation of internal-node values to be used as the server-domain
    3276 !--   values. This subroutine is based on the first-order numerical
    3277 !--   integration of the fine-grid values contained within the coarse-grid
    3278 !--   cell.
    3279 !     
    3280 !--                                        Antti Hellsten 16.9.2015.
    3281 !
    3282       IMPLICIT NONE
    3283       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN)  ::  f     !:
    3284       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)   ::  fc    !:
    3285       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifl   !:
    3286       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)                    ::  ifu   !:
    3287       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfl   !:
    3288       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)                    ::  jfu   !:
    3289       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)        ::  kb    !: may be unnecessary
    3290       INTEGER(iwp), INTENT(IN)                                        ::  kce   !:
    3291       INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfl   !:
    3292       INTEGER(iwp), DIMENSION(0:kce), INTENT(IN)                      ::  kfu   !:
    3293       CHARACTER(LEN=1), INTENT(IN)                                    ::  var   !:                           
    3294 
    3295       INTEGER(iwp) ::  i         !:
    3296       INTEGER(iwp) ::  icb       !:
    3297       INTEGER(iwp) ::  ice       !:
    3298       INTEGER(iwp) ::  ifc       !:
    3299       INTEGER(iwp) ::  ijfc      !:
    3300       INTEGER(iwp) ::  j         !:
    3301       INTEGER(iwp) ::  jcb       !:
    3302       INTEGER(iwp) ::  jce       !:
    3303       INTEGER(iwp) ::  k         !:
    3304       INTEGER(iwp) ::  kcb       !:
    3305       INTEGER(iwp) ::  l         !:
    3306       INTEGER(iwp) ::  m         !:
    3307       INTEGER(iwp) ::  n         !:
    3308       INTEGER(iwp) ::  nfc       !:
    3309       REAL(wp)     ::  cellsum   !:
    3310       REAL(wp)     ::  f1f       !:
    3311       REAL(wp)     ::  fra       !:
     3910       IF ( var == 'w' )  THEN
     3911          k    = nzt
     3912          ked  = nzt - 1
     3913       ELSE
     3914          k    = nzt + 1
     3915          ked  = nzt
     3916       ENDIF
     3917
     3918       DO  i = nxl, nxr
     3919          DO  j = nys, nyn
     3920             vdotnor = w(ked,j,i)
     3921!
     3922!--          Local outflow
     3923             IF ( vdotnor > 0.0_wp )  THEN
     3924                f(k,j,i) = f(ked,j,i)
     3925             ENDIF
     3926          ENDDO
     3927       ENDDO
     3928
     3929!
     3930!--    Just fill up the second ghost-node layer for w
     3931       IF ( var == 'w' )  THEN
     3932          f(nzt+1,:,:) = f(nzt,:,:)
     3933       ENDIF
     3934
     3935    END SUBROUTINE pmci_extrap_ifoutflow_t
     3936
     3937
     3938
     3939    SUBROUTINE pmci_anterp_tophat( f, fc, kct, ifl, ifu, jfl, jfu, kfl, kfu,   &
     3940                                   kb, var )
     3941!
     3942!--    Anterpolation of internal-node values to be used as the server-domain
     3943!--    values. This subroutine is based on the first-order numerical
     3944!--    integration of the fine-grid values contained within the coarse-grid
     3945!--    cell.
     3946       IMPLICIT NONE
     3947
     3948       CHARACTER(LEN=1), INTENT(IN) ::  var   !:
     3949
     3950       INTEGER(iwp) ::  i         !: Fine-grid index
     3951       INTEGER(iwp) ::  ii        !: Coarse-grid index
     3952       INTEGER(iwp) ::  iclp      !:
     3953       INTEGER(iwp) ::  icrm      !:
     3954       INTEGER(iwp) ::  ifc       !:
     3955       INTEGER(iwp) ::  ijfc      !:
     3956       INTEGER(iwp) ::  j         !: Fine-grid index
     3957       INTEGER(iwp) ::  jj        !: Coarse-grid index
     3958       INTEGER(iwp) ::  jcnm      !:
     3959       INTEGER(iwp) ::  jcsp      !:
     3960       INTEGER(iwp) ::  k         !: Fine-grid index
     3961       INTEGER(iwp) ::  kk        !: Coarse-grid index
     3962       INTEGER(iwp) ::  kcb       !:
     3963       INTEGER(iwp) ::  nfc       !:
     3964
     3965       INTEGER(iwp), INTENT(IN) ::  kct   !:
     3966
     3967       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)             ::  ifl   !:
     3968       INTEGER(iwp), DIMENSION(icl:icr), INTENT(IN)             ::  ifu   !:
     3969       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)             ::  jfl   !:
     3970       INTEGER(iwp), DIMENSION(jcs:jcn), INTENT(IN)             ::  jfu   !:
     3971!--    TO_DO: is the next line really unnecessary?
     3972       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN) ::  kb    !: may be unnecessary
     3973       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)               ::  kfl   !:
     3974       INTEGER(iwp), DIMENSION(0:kct), INTENT(IN)               ::  kfu   !:
     3975
     3976
     3977       REAL(wp) ::  cellsum   !:
     3978       REAL(wp) ::  f1f       !:
     3979       REAL(wp) ::  fra       !:
     3980
     3981       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  f   !:
     3982       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(INOUT)  ::  fc  !:
    33123983 
    33133984
    3314       icb = icl
    3315       ice = icr
    3316       jcb = jcs
    3317       jce = jcn
    3318 
    3319 !
    3320 !--   Define the index bounds icb, ice, jcb and jce.
    3321 !--   Note that kcb is simply zero and kce enters here as a parameter and it is
    3322 !--   determined in init_anterp_tophat 
    3323       IF ( nest_bound_l )  THEN
    3324          IF ( var == 'u' )  THEN
    3325             icb = icl + nhll + 1
    3326          ELSE
    3327             icb = icl + nhll
    3328          ENDIF
    3329       ENDIF
    3330       IF ( nest_bound_r )  THEN
    3331          ice = icr - nhlr
    3332       ENDIF
    3333 
    3334       IF ( nest_bound_s )  THEN
    3335          IF ( var == 'v' )  THEN
    3336             jcb = jcs + nhls + 1
    3337          ELSE
    3338             jcb = jcs + nhls
    3339          ENDIF
    3340       ENDIF
    3341       IF ( nest_bound_n )  THEN
    3342          jce = jcn - nhln
    3343       ENDIF
    3344       kcb = 0
    3345 
    3346 !
    3347 !--   Note that l,m, and n are coarse-grid indices and i,j, and k are fine-grid indices.
    3348       DO  l = icb, ice
    3349          ifc = ifu(l) - ifl(l) + 1
    3350          DO  m = jcb, jce
    3351             ijfc = ifc * ( jfu(m) - jfl(m) +1 )
    3352 
    3353 !
    3354 !--         How to deal with the lower bound of k in case of non-flat topography?
    3355             !kcb = MIN( kb(jfl(m),ifl(l)), kb(jfu(m),ifl(l)), kb(jfl(m),ifu(l)), kb(jfu(m),ifu(l)) ) ! Something wrong with this.
    3356             DO  n = kcb, kce
    3357                nfc =  ijfc * ( kfu(n) - kfl(n) + 1 )
    3358                cellsum = 0.0
    3359                DO  i = ifl(l), ifu(l)
    3360                   DO  j = jfl(m), jfu(m)
    3361                      DO  k = kfl(n), kfu(n)
    3362                         cellsum = cellsum + f(k,j,i)
    3363                      ENDDO
    3364                   ENDDO
    3365                ENDDO
    3366 
    3367 !
    3368 !--            Spatial under-relaxation.
    3369                fra  = frax(l) * fray(m) * fraz(n)
    3370 !--            TO_DO: why not KIND=wp ?
    3371                fc(n,m,l) = ( 1.0_wp - fra ) * fc(n,m,l) + fra * cellsum / REAL( nfc, KIND=KIND(cellsum) ) 
    3372             ENDDO
    3373          ENDDO
    3374       ENDDO
    3375 
    3376    END SUBROUTINE pmci_anterp_tophat
     3985!
     3986!--    Initialize the index bounds for anterpolation
     3987       iclp = icl
     3988       icrm = icr
     3989       jcsp = jcs
     3990       jcnm = jcn
     3991
     3992!
     3993!--    Define the index bounds iclp, icrm, jcsp and jcnm.
     3994!--    Note that kcb is simply zero and kct enters here as a parameter and it is
     3995!--    determined in pmci_init_anterp_tophat
     3996       IF ( nest_bound_l )  THEN
     3997          IF ( var == 'u' )  THEN
     3998             iclp = icl + nhll + 1
     3999          ELSE
     4000             iclp = icl + nhll
     4001          ENDIF
     4002       ENDIF
     4003       IF ( nest_bound_r )  THEN
     4004          icrm = icr - nhlr
     4005       ENDIF
     4006
     4007       IF ( nest_bound_s )  THEN
     4008          IF ( var == 'v' )  THEN
     4009             jcsp = jcs + nhls + 1
     4010          ELSE
     4011             jcsp = jcs + nhls
     4012          ENDIF
     4013       ENDIF
     4014       IF ( nest_bound_n )  THEN
     4015          jcnm = jcn - nhln
     4016       ENDIF
     4017       kcb = 0
     4018
     4019!
     4020!--    Note that l,m, and n are coarse-grid indices and i,j, and k are fine-grid
     4021!--    indices.
     4022       DO  ii = iclp, icrm
     4023          ifc = ifu(ii) - ifl(ii) + 1
     4024          DO  jj = jcsp, jcnm
     4025             ijfc = ifc * ( jfu(jj) - jfl(jj) + 1 )
     4026!
     4027!--          For simplicity anterpolate within buildings too
     4028             DO  kk = kcb, kct
     4029                nfc =  ijfc * ( kfu(kk) - kfl(kk) + 1 )
     4030                cellsum = 0.0_wp
     4031                DO  i = ifl(ii), ifu(ii)
     4032                   DO  j = jfl(jj), jfu(jj)
     4033                      DO  k = kfl(kk), kfu(kk)
     4034                         cellsum = cellsum + f(k,j,i)
     4035                      ENDDO
     4036                   ENDDO
     4037                ENDDO
     4038!
     4039!--             Spatial under-relaxation.
     4040                fra  = frax(ii) * fray(jj) * fraz(kk)
     4041!
     4042!--             Block out the fine-grid corner patches from the anterpolation
     4043                IF ( ( ifl(ii) < nxl ) .OR. ( ifu(ii) > nxr ) )  THEN
     4044                   IF ( ( jfl(jj) < nys ) .OR. ( jfu(jj) > nyn ) )  THEN
     4045                      fra = 0.0_wp
     4046                   ENDIF
     4047                ENDIF
     4048
     4049                fc(kk,jj,ii) = ( 1.0_wp - fra ) * fc(kk,jj,ii) +                 &
     4050                               fra * cellsum / REAL( nfc, KIND = wp )
     4051
     4052             ENDDO
     4053          ENDDO
     4054       ENDDO
     4055
     4056    END SUBROUTINE pmci_anterp_tophat
    33774057
    33784058#endif
    33794059 END SUBROUTINE pmci_client_datatrans
    33804060
    3381 
    3382 
    3383  SUBROUTINE pmci_update_new
    3384 
    3385 #if defined( __parallel )
    3386 !
    3387 !-- Copy the interpolated/anterpolated boundary values to the _p
    3388 !-- arrays, too, to make sure the interpolated/anterpolated boundary
    3389 !-- values are carried over from one RK inner step to another.
    3390 !-- So far works only with the cpp-switch __nopointer.
    3391 !
    3392 !--                 Antti Hellsten 8.3.2015
    3393 !
    3394 
    3395 !-- Just debugging
    3396     w(nzt+1,:,:) = w(nzt,:,:)
    3397 
    3398     u_p  = u
    3399     v_p  = v
    3400     w_p  = w
    3401     e_p  = e
    3402     pt_p = pt
    3403     IF ( humidity  .OR.  passive_scalar )  THEN
    3404        q_p = q
    3405     ENDIF
    3406 
    3407 !
    3408 !-- TO_DO: Find out later if nesting would work without __nopointer.
    3409 #endif
    3410 
    3411  END SUBROUTINE pmci_update_new
    3412 
    3413 
    3414 
    3415  SUBROUTINE pmci_set_array_pointer( name, client_id, nz_cl )
    3416 
    3417     IMPLICIT NONE
    3418 
    3419     INTEGER, INTENT(IN)                  ::  client_id   !:
    3420     INTEGER, INTENT(IN)                  ::  nz_cl       !:
    3421     CHARACTER(LEN=*), INTENT(IN)         ::  name        !:
    3422    
    3423 #if defined( __parallel )
    3424     REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d        !:
    3425     REAL(wp), POINTER, DIMENSION(:,:)    ::  p_2d_sec    !:
    3426     REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d        !:
    3427     REAL(wp), POINTER, DIMENSION(:,:,:)  ::  p_3d_sec    !:
    3428     INTEGER(iwp)                         ::  ierr        !:
    3429     INTEGER(iwp)                         ::  istat       !:
    3430    
    3431 
    3432     NULLIFY( p_3d )
    3433     NULLIFY( p_2d )
    3434 
    3435 !
    3436 !-- List of array names, which can be coupled
    3437 !-- In case of 3D please change also the second array for the pointer version
    3438     IF ( TRIM(name) == "u" )     p_3d => u
    3439     IF ( TRIM(name) == "v" )     p_3d => v
    3440     IF ( TRIM(name) == "w" )     p_3d => w
    3441     IF ( TRIM(name) == "e" )     p_3d => e
    3442     IF ( TRIM(name) == "pt" )    p_3d => pt
    3443     IF ( TRIM(name) == "q" )     p_3d => q
    3444 !
    3445 !-- This is just an example for a 2D array, not active for coupling
    3446 !-- Please note, that z0 has to be declared as TARGET array in modules.f90
    3447 !    IF ( TRIM(name) == "z0" )    p_2d => z0
    3448 
    3449 #if defined( __nopointer )
    3450     IF ( ASSOCIATED( p_3d ) )  THEN
    3451        CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz )
    3452     ELSEIF ( ASSOCIATED( p_2d ) )  THEN
    3453        CALL pmc_s_set_dataarray( client_id, p_2d )
    3454     ELSE
    3455 !
    3456 !--    Give only one message for the root domain
    3457        IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
    3458 
    3459           message_string = 'pointer for array "' // TRIM( name ) //            &
    3460                            '" can''t be associated'
    3461           CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
    3462        ELSE
    3463 !
    3464 !--       Avoid others to continue
    3465           CALL MPI_BARRIER( comm2d, ierr )
    3466        ENDIF
    3467     ENDIF
    3468 #else
    3469 !-- TO_DO: Why aren't the other pointers (p_3d) not set to u_1, v_1, etc.??
    3470     IF ( TRIM(name) == "u" )     p_3d_sec => u_2
    3471     IF ( TRIM(name) == "v" )     p_3d_sec => v_2
    3472     IF ( TRIM(name) == "w" )     p_3d_sec => w_2
    3473     IF ( TRIM(name) == "e" )     p_3d_sec => e_2
    3474     IF ( TRIM(name) == "pt" )    p_3d_sec => pt_2
    3475     IF ( TRIM(name) == "q" )     p_3d_sec => q_2
    3476 
    3477     IF ( ASSOCIATED( p_3d ) )  THEN
    3478        CALL pmc_s_set_dataarray( client_id, p_3d, nz_cl, nz, &
    3479                                  array_2 = p_3d_sec )
    3480     ELSEIF ( ASSOCIATED( p_2d ) )  THEN
    3481        CALL pmc_s_set_dataarray( client_id, p_2d )
    3482     ELSE
    3483 !
    3484 !--    Give only one message for the root domain
    3485        IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
    3486 
    3487           message_string = 'pointer for array "' // TRIM( name ) //            &
    3488                            '" can''t be associated'
    3489           CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
    3490        ELSE
    3491 !
    3492 !--       Avoid others to continue
    3493           CALL MPI_BARRIER( comm2d, ierr )
    3494        ENDIF
    3495 
    3496    ENDIF
    3497 #endif
    3498 
    3499 #endif
    3500  END SUBROUTINE pmci_set_array_pointer
    3501 
    3502 
    3503 
    3504  SUBROUTINE pmci_create_client_arrays( name, is, ie, js, je, nzc  )
    3505 
    3506     IMPLICIT NONE
    3507 
    3508     INTEGER(iwp), INTENT(IN)              ::  ie      !:
    3509     INTEGER(iwp), INTENT(IN)              ::  is      !:
    3510     INTEGER(iwp), INTENT(IN)              ::  je      !:
    3511     INTEGER(iwp), INTENT(IN)              ::  js      !:
    3512     INTEGER(iwp), INTENT(IN)              ::  nzc     !:  Note that nzc is cg%nz
    3513     CHARACTER(LEN=*), INTENT(IN)          ::  name    !:
    3514      
    3515 #if defined( __parallel )
    3516     REAL(wp), POINTER,DIMENSION(:,:)      ::  p_2d    !:
    3517     REAL(wp), POINTER,DIMENSION(:,:,:)    ::  p_3d    !:
    3518     INTEGER(iwp)                          ::  ierr    !:
    3519     INTEGER(iwp)                          ::  istat   !:
    3520 
    3521 
    3522     NULLIFY( p_3d )
    3523     NULLIFY( p_2d )
    3524 
    3525 !
    3526 !-- List of array names, which can be coupled.
    3527 !-- AH: Note that the k-range of the *c arrays is changed from 1:nz to 0:nz+1.
    3528     IF ( TRIM( name ) == "u" )  THEN
    3529        IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1, js:je, is:ie) )
    3530        p_3d => uc
    3531     ELSEIF ( TRIM( name ) == "v" )  THEN
    3532        IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1, js:je, is:ie) )
    3533        p_3d => vc
    3534     ELSEIF ( TRIM( name ) == "w" )  THEN
    3535        IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1, js:je, is:ie) )
    3536        p_3d => wc
    3537     ELSEIF ( TRIM( name ) == "e" )  THEN
    3538        IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1, js:je, is:ie) )
    3539        p_3d => ec
    3540     ELSEIF ( TRIM( name ) == "pt")  THEN
    3541        IF ( .NOT. ALLOCATED( ptc ) ) ALLOCATE( ptc(0:nzc+1, js:je, is:ie) )
    3542        p_3d => ptc
    3543     ELSEIF ( TRIM( name ) == "q")  THEN
    3544        IF ( .NOT. ALLOCATED( qc ) ) ALLOCATE( qc(0:nzc+1, js:je, is:ie) )
    3545        p_3d => qc
    3546     !ELSEIF (trim(name) == "z0") then
    3547        !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
    3548        !p_2d => z0c
    3549     ENDIF
    3550 
    3551     IF ( ASSOCIATED( p_3d ) )  THEN
    3552        CALL pmc_c_set_dataarray( p_3d )
    3553     ELSEIF ( ASSOCIATED( p_2d ) )  THEN
    3554        CALL pmc_c_set_dataarray( p_2d )
    3555     ELSE
    3556 !
    3557 !--    Give only one message for the first client domain
    3558        IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
    3559 
    3560           message_string = 'pointer for array "' // TRIM( name ) //            &
    3561                            '" can''t be associated'
    3562           CALL message( 'pmci_create_client_arrays', 'PA0170', 3, 2, 0, 6, 0 )
    3563        ELSE
    3564 !
    3565 !--       Avoid others to continue
    3566           CALL MPI_BARRIER( comm2d, ierr )
    3567        ENDIF
    3568     ENDIF
    3569 
    3570 #endif
    3571  END SUBROUTINE pmci_create_client_arrays
    3572 
    3573 
    3574 
    3575  SUBROUTINE pmci_server_initialize
    3576 
    3577 #if defined( __parallel )
    3578     IMPLICIT NONE
    3579 
    3580     INTEGER(iwp)   ::  client_id   !:
    3581     INTEGER(iwp)   ::  m           !:
    3582     REAL(wp)       ::  waittime    !:
    3583 
    3584 
    3585     DO  m = 1, SIZE( pmc_server_for_client ) - 1
    3586        client_id = pmc_server_for_client(m)
    3587        CALL pmc_s_fillbuffer( client_id, waittime=waittime )
    3588     ENDDO
    3589 
    3590 #endif
    3591  END SUBROUTINE pmci_server_initialize
    3592 
    3593 
    3594 
    3595  SUBROUTINE pmci_client_initialize
    3596 
    3597 #if defined( __parallel )
    3598     IMPLICIT NONE
    3599 
    3600     INTEGER(iwp)   ::  i          !:
    3601     INTEGER(iwp)   ::  icl        !:
    3602     INTEGER(iwp)   ::  icr        !:
    3603     INTEGER(iwp)   ::  j          !:
    3604     INTEGER(iwp)   ::  jcn        !:
    3605     INTEGER(iwp)   ::  jcs        !:
    3606     REAL(wp)       ::  waittime   !:
    3607 
    3608 
    3609     IF ( cpl_id > 1 )  THEN   !  Root id is never a client     
    3610 
    3611 !
    3612 !--    Client domain boundaries in the server indice space.
    3613        icl = coarse_bound(1)
    3614        icr = coarse_bound(2)
    3615        jcs = coarse_bound(3)
    3616        jcn = coarse_bound(4)
    3617 
    3618 !
    3619 !--    Get data from server     
    3620        CALL pmc_c_getbuffer( waittime = waittime )
    3621 
    3622 !
    3623 !--    The interpolation.
    3624        CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo, r2yo, r1zo, r2zo, nzb_u_inner, 'u' )
    3625        CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv, r2yv, r1zo, r2zo, nzb_v_inner, 'v' )
    3626        CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo, r2yo, r1zw, r2zw, nzb_w_inner, 'w' )
    3627        CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 'e' )
    3628        CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
    3629        IF ( humidity  .OR.  passive_scalar )  THEN
    3630           CALL pmci_interp_tril_all ( q, qc, ico, jco, kco, r1xo, r2xo, r1yo, r2yo, r1zo, r2zo, nzb_s_inner, 's' )
    3631        ENDIF
    3632 
    3633        IF ( topography /= 'flat' )  THEN
    3634 
    3635 !
    3636 !--       Inside buildings set velocities and TKE back to zero.
    3637 !--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
    3638 !--       maybe revise later.
    3639           DO   i = nxlg, nxrg
    3640              DO   j = nysg, nyng
    3641                 u(nzb:nzb_u_inner(j,i),j,i)   = 0.0_wp
    3642                 v(nzb:nzb_v_inner(j,i),j,i)   = 0.0_wp
    3643                 w(nzb:nzb_w_inner(j,i),j,i)   = 0.0_wp
    3644                 e(nzb:nzb_s_inner(j,i),j,i)   = 0.0_wp
    3645                 u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0_wp
    3646                 v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0_wp
    3647                 w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0_wp
    3648                 e_p(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
    3649              ENDDO
    3650           ENDDO
    3651        ENDIF
    3652     ENDIF
    3653    
    3654 
    3655  CONTAINS
    3656 
    3657 
    3658    SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, r2z, kb, var )
    3659 
    3660 !
    3661 !-- Interpolation of the internal values for the client-domain initialization.
    3662 !-- This subroutine is based on trilinear interpolation.
    3663 !-- Coding based on interp_tril_lr/sn/t
    3664 !
    3665 !--                                     Antti Hellsten 20.10.2015.
    3666       IMPLICIT NONE
    3667       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) ::  f     !:
    3668       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN)        ::  fc    !:
    3669       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r1x   !:
    3670       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)                        ::  r2x   !:
    3671       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r1y   !:
    3672       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)                        ::  r2y   !:
    3673       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r1z   !:
    3674       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)                        ::  r2z   !:
    3675       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)                    ::  ic    !:
    3676       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)                    ::  jc    !:
    3677       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)                    ::  kc    !:
    3678       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg), INTENT(IN)          ::  kb    !:
    3679       CHARACTER(LEN=1), INTENT(IN) :: var  !:
    3680      
    3681       INTEGER(iwp) ::  i      !:
    3682       INTEGER(iwp) ::  ib     !:
    3683       INTEGER(iwp) ::  ie     !:
    3684       INTEGER(iwp) ::  j      !:
    3685       INTEGER(iwp) ::  jb     !:
    3686       INTEGER(iwp) ::  je     !:
    3687       INTEGER(iwp) ::  k      !:
    3688       INTEGER(iwp) ::  k1     !:
    3689       INTEGER(iwp) ::  kbc    !:
    3690       INTEGER(iwp) ::  l      !:
    3691       INTEGER(iwp) ::  m      !:
    3692       INTEGER(iwp) ::  n      !:
    3693       REAL(wp) ::  fk         !:
    3694       REAL(wp) ::  fkj        !:
    3695       REAL(wp) ::  fkjp       !:
    3696       REAL(wp) ::  fkp        !:
    3697       REAL(wp) ::  fkpj       !:
    3698       REAL(wp) ::  fkpjp      !:
    3699       REAL(wp) ::  logratio   !:
    3700       REAL(wp) ::  logzuc1    !:
    3701       REAL(wp) ::  zuc1       !:
    3702      
    3703      
    3704       ib = nxl
    3705       ie = nxr
    3706       jb = nys   
    3707       je = nyn
    3708       IF ( nest_bound_l )  THEN
    3709          ib = nxl - 1
    3710          IF ( var == 'u' )  THEN   !  For u, nxl is a ghost node, but not for the other variables.
    3711             ib = nxl
    3712          ENDIF
    3713       ENDIF
    3714       IF ( nest_bound_s )  THEN
    3715          jb = nys - 1
    3716          IF ( var == 'v' )  THEN   !  For v, nys is a ghost node, but not for the other variables.
    3717             jb = nys
    3718          ENDIF
    3719       ENDIF
    3720       IF ( nest_bound_r )  THEN
    3721          ie = nxr + 1
    3722       ENDIF
    3723       IF ( nest_bound_n )  THEN
    3724          je = nyn + 1
    3725       ENDIF
    3726 
    3727 !
    3728 !--   Trilinear interpolation.
    3729       DO  i = ib, ie
    3730          DO  j = jb, je
    3731             DO  k = kb(j,i), nzt + 1
    3732                l = ic(i)
    3733                m = jc(j)
    3734                n = kc(k)
    3735                fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
    3736                fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
    3737                fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
    3738                fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
    3739                fk       = r1y(j) * fkj  + r2y(j) * fkjp
    3740                fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
    3741                f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
    3742             ENDDO
    3743          ENDDO
    3744       ENDDO
    3745 
    3746 !
    3747 !--   Correct the interpolated values of u and v in near-wall nodes, i.e. in
    3748 !--   the nodes below the coarse-grid nodes with k=1. The corrction is only made
    3749 !--   over horizontal wall surfaces in this phase. For the nest boundary conditions,
    3750 !--   a corresponding corrections is made for all vertical walls, too.
    3751       IF ( var == 'u' .OR. var == 'v' )  THEN
    3752          DO  i = ib, nxr
    3753             DO  j = jb, nyn             
    3754                kbc = 1
    3755                DO  WHILE ( cg%zu(kbc) < zu(kb(j,i)) )   !  kbc is the first coarse-grid point above the surface.
    3756                   kbc = kbc + 1
    3757                ENDDO
    3758                zuc1    = cg%zu(kbc)
    3759                k1      = kb(j,i) + 1 
    3760                DO  WHILE ( zu(k1) < zuc1 )
    3761                   k1 = k1 + 1
    3762                ENDDO
    3763                logzuc1 = LOG( ( zu(k1) - zu(kb(j,i)) ) / z0(j,i) )
    3764                
    3765                k = kb(j,i) + 1
    3766                DO  WHILE ( zu(k) < zuc1 )
    3767                   logratio = ( LOG( ( zu(k) - zu(kb(j,i)) ) / z0(j,i)) ) / logzuc1
    3768                   f(k,j,i) = logratio * f(k1,j,i)
    3769                   k  = k + 1
    3770                ENDDO
    3771                f(kb(j,i),j,i) = 0.0_wp
    3772             ENDDO
    3773          ENDDO
    3774       ELSEIF ( var == 'w' )  THEN
    3775          DO  i = ib, nxr
    3776             DO  j = jb, nyn
    3777                f(kb(j,i),j,i) = 0.0_wp
    3778             ENDDO
    3779          ENDDO
    3780       ENDIF
    3781      
    3782    END SUBROUTINE pmci_interp_tril_all
    3783 
    3784 #endif
    3785  END SUBROUTINE pmci_client_initialize
    3786 
    3787 
    3788 
    3789  SUBROUTINE pmci_ensure_nest_mass_conservation
    3790 
    3791 #if defined( __parallel )
    3792 !
    3793 !-- Adjust the volume-flow rate through the top boundary
    3794 !-- so that the net volume flow through all boundaries
    3795 !-- of the current nest domain becomes zero.
    3796     IMPLICIT NONE
    3797 
    3798     INTEGER(iwp)  ::  i                          !:
    3799     INTEGER(iwp)  ::  ierr                       !:
    3800     INTEGER(iwp)  ::  j                          !:
    3801     INTEGER(iwp)  ::  k                          !:
    3802     REAL(wp) ::  dxdy                            !:
    3803     REAL(wp) ::  innor                           !:
    3804     REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !:
    3805     REAL(wp) ::  w_lt                            !:
    3806 
    3807 !
    3808 !-- Sum up the volume flow through the left/right boundaries.
    3809     volume_flow(1)   = 0.0_wp
    3810     volume_flow_l(1) = 0.0_wp
    3811 
    3812     IF ( nest_bound_l )  THEN
    3813        i = 0
    3814        innor = dy
    3815        DO   j = nys, nyn
    3816           DO   k = nzb_u_inner(j,i) + 1, nzt
    3817              volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)           
    3818           ENDDO
    3819        ENDDO
    3820     ENDIF
    3821 
    3822     IF ( nest_bound_r )  THEN     
    3823        i = nx + 1
    3824        innor = -dy
    3825        DO   j = nys, nyn
    3826           DO   k = nzb_u_inner(j,i) + 1, nzt
    3827              volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)           
    3828           ENDDO
    3829        ENDDO
    3830     ENDIF
    3831      
    3832 #if defined( __parallel )   
    3833     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    3834     CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
    3835                         MPI_SUM, comm2d, ierr )
    3836 #else
    3837     volume_flow(1) = volume_flow_l(1)
    3838 #endif
    3839 
    3840 !
    3841 !-- Sum up the volume flow through the south/north boundaries.
    3842     volume_flow(2)   = 0.0_wp
    3843     volume_flow_l(2) = 0.0_wp
    3844 
    3845     IF ( nest_bound_s )  THEN
    3846        j = 0
    3847        innor = dx
    3848        DO   i = nxl, nxr
    3849           DO   k = nzb_v_inner(j,i) + 1, nzt
    3850              volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
    3851           ENDDO
    3852        ENDDO
    3853     ENDIF
    3854 
    3855     IF ( nest_bound_n )  THEN
    3856        j = ny + 1
    3857        innor = -dx
    3858        DO   i = nxl, nxr
    3859           DO   k = nzb_v_inner(j,i)+1, nzt
    3860              volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)
    3861           ENDDO
    3862        ENDDO
    3863     ENDIF
    3864 
    3865 #if defined( __parallel )   
    3866     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    3867     CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
    3868                         MPI_SUM, comm2d, ierr )   
    3869 #else
    3870     volume_flow(2) = volume_flow_l(2)
    3871 #endif
    3872      
    3873 !
    3874 !-- Sum up the volume flow through the top boundary.
    3875     volume_flow(3)   = 0.0_wp
    3876     volume_flow_l(3) = 0.0_wp
    3877     dxdy = dx * dy
    3878     k = nzt
    3879     DO   i = nxl, nxr
    3880        DO   j = nys, nyn
    3881           volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
    3882        ENDDO
    3883     ENDDO
    3884        
    3885 #if defined( __parallel )   
    3886     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    3887     CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL, &
    3888                         MPI_SUM, comm2d, ierr )   
    3889 #else
    3890     volume_flow(3) = volume_flow_l(3)
    3891 #endif
    3892 
    3893 !
    3894 !-- Correct the top-boundary value of w.
    3895     w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
    3896     DO   i = nxl, nxr
    3897        DO   j = nys, nyn
    3898           DO  k = nzt, nzt + 1
    3899              w(k,j,i) = w(k,j,i) + w_lt
    3900           ENDDO
    3901        ENDDO
    3902     ENDDO
    3903 
    3904 #endif
    3905  END SUBROUTINE pmci_ensure_nest_mass_conservation
    3906 
    39074061END MODULE pmc_interface
  • palm/trunk/SOURCE/pmc_server.f90

    r1787 r1791  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Debug write-statements commented out
    2323!
    2424! Former revisions:
     
    136136
    137137      do i=1,size(PMC_Server_for_Client)-1
    138          if(m_model_comm == 0) write(0,*) 'PMC_Server: Initialize client Id',PMC_Server_for_Client(i)
     138!         if(m_model_comm == 0) write(0,*) 'PMC_Server: Initialize client Id',PMC_Server_for_Client(i)
    139139
    140140         ClientId = PMC_Server_for_Client(i)
     
    152152         CALL MPI_Comm_rank (Clients(ClientId)%intra_comm, Clients(ClientId)%intra_rank, istat);
    153153
    154          write(9,*) 'ClientId ',i,ClientId,m_world_rank, Clients(ClientId)%inter_npes
     154!         write(9,*) 'ClientId ',i,ClientId,m_world_rank, Clients(ClientId)%inter_npes
    155155
    156156         ALLOCATE (Clients(ClientId)%PEs(Clients(ClientId)%inter_npes))
  • palm/trunk/SOURCE/time_integration.f90

    r1787 r1791  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! call of pmci_update_new removed
    2222!
    2323! Former revisions:
     
    257257               pmci_ensure_nest_mass_conservation, pmci_client_datatrans,      &
    258258               pmci_client_synchronize, pmci_server_datatrans,                 &
    259                pmci_server_synchronize, pmci_update_new, server_to_client
     259               pmci_server_synchronize, server_to_client
    260260
    261261    USE production_e_mod,                                                      &
     
    703703             IF ( nest_domain )  THEN
    704704                CALL pmci_ensure_nest_mass_conservation
    705 !
    706 !--             pmc_update_new is not necessary if nesting is made at each
    707 !--             substep
    708                 CALL pmci_update_new
    709705             ENDIF
    710706
Note: See TracChangeset for help on using the changeset viewer.