source: palm/trunk/SOURCE/pmc_interface_mod.f90 @ 2951

Last change on this file since 2951 was 2951, checked in by kanani, 3 years ago

Add cpu measures for nesting initialization

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/forwind/SOURCE/pmc_interface_mod.f901564-1913
    /palm/branches/palm4u/SOURCE/pmc_interface_mod.f902540-2692
    /palm/trunk/SOURCE/pmc_interface_mod.f90mergedeligible
    /palm/branches/fricke/SOURCE/pmc_interface_mod.f90942-977
    /palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90989-1052
    /palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90296-409
    /palm/branches/suehring/pmc_interface_mod.f90423-666
File size: 221.6 KB
Line 
1!> @file pmc_interface_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: pmc_interface_mod.f90 2951 2018-04-06 09:05:08Z kanani $
27! Add log_point_s for pmci_model_configuration
28!
29! 2938 2018-03-27 15:52:42Z suehring
30! - Nesting for RANS mode implemented
31!    - Interpolation of TKE onto child domain only if both parent and child are
32!      either in LES mode or in RANS mode
33!    - Interpolation of dissipation if both parent and child are in RANS mode
34!      and TKE-epsilon closure is applied
35!    - Enable anterpolation of TKE and dissipation rate in case parent and
36!      child operate in RANS mode
37!
38! - Some unused variables removed from ONLY list
39! - Some formatting adjustments for particle nesting
40!
41! 2936 2018-03-27 14:49:27Z suehring
42! Control logics improved to allow nesting also in cases with
43! constant_flux_layer = .F. or constant_diffusion = .T.
44!
45! 2895 2018-03-15 10:26:12Z hellstea
46! Change in the nest initialization (pmci_interp_tril_all). Bottom wall BC is no
47! longer overwritten.
48!
49! 2868 2018-03-09 13:25:09Z hellstea
50! Local conditional Neumann conditions for one-way coupling removed. 
51!
52! 2853 2018-03-05 14:44:20Z suehring
53! Bugfix in init log-law correction.
54!
55! 2841 2018-02-27 15:02:57Z knoop
56! Bugfix: wrong placement of include 'mpif.h' corrected
57!
58! 2812 2018-02-16 13:40:25Z hellstea
59! Bugfixes in computation of the interpolation loglaw-correction parameters
60!
61! 2809 2018-02-15 09:55:58Z schwenkel
62! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
63!
64! 2806 2018-02-14 17:10:37Z thiele
65! Bugfixing Write statements
66!
67! 2804 2018-02-14 16:57:03Z thiele
68! preprocessor directive for c_sizeof in case of __gfortran added
69!
70! 2803 2018-02-14 16:56:32Z thiele
71! Introduce particle transfer in nested models.
72!
73! 2795 2018-02-07 14:48:48Z hellstea
74! Bugfix in computation of the anterpolation under-relaxation functions.
75!
76! 2773 2018-01-30 14:12:54Z suehring
77! - Nesting for chemical species
78! - Bugfix in setting boundary condition at downward-facing walls for passive
79!   scalar
80! - Some formatting adjustments
81!
82! 2718 2018-01-02 08:49:38Z maronga
83! Corrected "Former revisions" section
84!
85! 2701 2017-12-15 15:40:50Z suehring
86! Changes from last commit documented
87!
88! 2698 2017-12-14 18:46:24Z suehring
89! Bugfix in get_topography_top_index
90!
91! 2696 2017-12-14 17:12:51Z kanani
92! Change in file header (GPL part)
93! - Bugfix in init_tke_factor (MS)
94!
95! 2669 2017-12-06 16:03:27Z raasch
96! file extension for nested domains changed to "_N##",
97! created flag file in order to enable combine_plot_fields to process nest data
98!
99! 2663 2017-12-04 17:40:50Z suehring
100! Bugfix, wrong coarse-grid index in init_tkefactor used.
101!
102! 2602 2017-11-03 11:06:41Z hellstea
103! Index-limit related bug (occurred with nesting_mode='vertical') fixed in
104! pmci_interp_tril_t. Check for too high nest domain added in pmci_setup_parent.   
105! Some cleaning up made.
106!
107! 2582 2017-10-26 13:19:46Z hellstea
108! Resetting of e within buildings / topography in pmci_parent_datatrans removed
109! as unnecessary since e is not anterpolated, and as incorrect since it overran
110! the default Neumann condition (bc_e_b).
111!
112! 2359 2017-08-21 07:50:30Z hellstea
113! A minor indexing error in pmci_init_loglaw_correction is corrected.
114!
115! 2351 2017-08-15 12:03:06Z kanani
116! Removed check (PA0420) for nopointer version
117!
118! 2350 2017-08-15 11:48:26Z kanani
119! Bugfix and error message for nopointer version.
120!
121! 2318 2017-07-20 17:27:44Z suehring
122! Get topography top index via Function call
123!
124! 2317 2017-07-20 17:27:19Z suehring
125! Set bottom boundary condition after anterpolation.
126! Some variable description added.
127!
128! 2293 2017-06-22 12:59:12Z suehring
129! In anterpolation, exclude grid points which are used for interpolation.
130! This avoids the accumulation of numerical errors leading to increased
131! variances for shallow child domains. 
132!
133! 2292 2017-06-20 09:51:42Z schwenkel
134! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
135! includes two more prognostic equations for cloud drop concentration (nc) 
136! and cloud water content (qc).
137!
138! 2285 2017-06-15 13:15:41Z suehring
139! Consider multiple pure-vertical nest domains in overlap check
140!
141! 2283 2017-06-14 10:17:34Z suehring
142! Bugfixes in inititalization of log-law correction concerning
143! new topography concept
144!
145! 2281 2017-06-13 11:34:50Z suehring
146! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
147!
148! 2241 2017-06-01 13:46:13Z hellstea
149! A minor indexing error in pmci_init_loglaw_correction is corrected.
150!
151! 2240 2017-06-01 13:45:34Z hellstea
152!
153! 2232 2017-05-30 17:47:52Z suehring
154! Adjustments to new topography concept
155!
156! 2229 2017-05-30 14:52:52Z hellstea
157! A minor indexing error in pmci_init_anterp_tophat is corrected.
158!
159! 2174 2017-03-13 08:18:57Z maronga
160! Added support for cloud physics quantities, syntax layout improvements. Data
161! transfer of qc and nc is prepared but currently deactivated until both
162! quantities become prognostic variables.
163! Some bugfixes.
164!
165! 2019 2016-09-30 13:40:09Z hellstea
166! Bugfixes mainly in determining the anterpolation index bounds. These errors
167! were detected when first time tested using 3:1 grid-spacing.
168!
169! 2003 2016-08-24 10:22:32Z suehring
170! Humidity and passive scalar also separated in nesting mode
171!
172! 2000 2016-08-20 18:09:15Z knoop
173! Forced header and separation lines into 80 columns
174!
175! 1938 2016-06-13 15:26:05Z hellstea
176! Minor clean-up.
177!
178! 1901 2016-05-04 15:39:38Z raasch
179! Initial version of purely vertical nesting introduced.
180! Code clean up. The words server/client changed to parent/child.
181!
182! 1900 2016-05-04 15:27:53Z raasch
183! unused variables removed
184!
185! 1894 2016-04-27 09:01:48Z raasch
186! bugfix: pt interpolations are omitted in case that the temperature equation is
187! switched off
188!
189! 1892 2016-04-26 13:49:47Z raasch
190! bugfix: pt is not set as a data array in case that the temperature equation is
191! switched off with neutral = .TRUE.
192!
193! 1882 2016-04-20 15:24:46Z hellstea
194! The factor ijfc for nfc used in anterpolation is redefined as 2-D array
195! and precomputed in pmci_init_anterp_tophat.
196!
197! 1878 2016-04-19 12:30:36Z hellstea
198! Synchronization rewritten, logc-array index order changed for cache
199! optimization
200!
201! 1850 2016-04-08 13:29:27Z maronga
202! Module renamed
203!
204!
205! 1808 2016-04-05 19:44:00Z raasch
206! MPI module used by default on all machines
207!
208! 1801 2016-04-05 13:12:47Z raasch
209! bugfix for r1797: zero setting of temperature within topography does not work
210! and has been disabled
211!
212! 1797 2016-03-21 16:50:28Z raasch
213! introduction of different datatransfer modes,
214! further formatting cleanup, parameter checks added (including mismatches
215! between root and nest model settings),
216! +routine pmci_check_setting_mismatches
217! comm_world_nesting introduced
218!
219! 1791 2016-03-11 10:41:25Z raasch
220! routine pmci_update_new removed,
221! pmc_get_local_model_info renamed pmc_get_model_info, some keywords also
222! renamed,
223! filling up redundant ghost points introduced,
224! some index bound variables renamed,
225! further formatting cleanup
226!
227! 1783 2016-03-06 18:36:17Z raasch
228! calculation of nest top area simplified,
229! interpolation and anterpolation moved to seperate wrapper subroutines
230!
231! 1781 2016-03-03 15:12:23Z raasch
232! _p arrays are set zero within buildings too, t.._m arrays and respective
233! settings within buildings completely removed
234!
235! 1779 2016-03-03 08:01:28Z raasch
236! only the total number of PEs is given for the domains, npe_x and npe_y
237! replaced by npe_total, two unused elements removed from array
238! parent_grid_info_real,
239! array management changed from linked list to sequential loop
240!
241! 1766 2016-02-29 08:37:15Z raasch
242! modifications to allow for using PALM's pointer version,
243! +new routine pmci_set_swaplevel
244!
245! 1764 2016-02-28 12:45:19Z raasch
246! +cpl_parent_id,
247! cpp-statements for nesting replaced by __parallel statements,
248! errors output with message-subroutine,
249! index bugfixes in pmci_interp_tril_all,
250! some adjustments to PALM style
251!
252! 1762 2016-02-25 12:31:13Z hellstea
253! Initial revision by A. Hellsten
254!
255! Description:
256! ------------
257! Domain nesting interface routines. The low-level inter-domain communication   
258! is conducted by the PMC-library routines.
259!
260! @todo Remove array_3d variables from USE statements thate not used in the
261!       routine
262! @todo Data transfer of qc and nc is prepared but not activated
263!------------------------------------------------------------------------------!
264 MODULE pmc_interface
265
266    USE ISO_C_BINDING
267
268
269#if defined( __nopointer )
270    USE arrays_3d,                                                             &
271        ONLY:  diss, dzu, dzw, e, e_p, nc, nr, pt, q, qc, qr, s, u, u_p,       &
272               v, v_p, w, w_p, zu, zw
273#else
274   USE arrays_3d,                                                              &
275        ONLY:  diss, diss_2, dzu, dzw, e, e_p, e_2, nc, nc_2, nc_p, nr, nr_2,  &
276               pt, pt_2, q, q_2, qc, qc_2, qr, qr_2, s, s_2,                       &
277               u, u_p, u_2, v, v_p, v_2, w, w_p, w_2, zu, zw
278#endif
279
280    USE control_parameters,                                                    &
281        ONLY:  air_chemistry, cloud_physics,                                   &
282               constant_diffusion, constant_flux_layer,                        &
283               coupling_char, dt_3d, dz, humidity, message_string,             &
284               microphysics_morrison, microphysics_seifert,                    &
285               nest_bound_l, nest_bound_r, nest_bound_s, nest_bound_n,         &
286               nest_domain, neutral, passive_scalar, rans_mode, rans_tke_e,    &
287               roughness_length, simulated_time, topography, volume_flow
288
289    USE chem_modules,                                                          &
290        ONLY:  nspec
291
292    USE chemistry_model_mod,                                                   &
293        ONLY:  chem_species, spec_conc_2
294
295    USE cpulog,                                                                &
296        ONLY:  cpu_log, log_point_s
297
298    USE grid_variables,                                                        &
299        ONLY:  dx, dy
300
301    USE indices,                                                               &
302        ONLY:  nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
303               nysv, nz, nzb, nzt, wall_flags_0
304
305    USE particle_attributes,                                                   &
306        ONLY:  particle_advection
307
308    USE kinds
309
310#if defined( __parallel )
311#if !defined( __mpifh )
312    USE MPI
313#endif
314
315    USE pegrid,                                                                &
316        ONLY:  collective_wait, comm1dx, comm1dy, comm2d, myid, myidx, myidy,  &
317               numprocs
318
319    USE pmc_child,                                                             &
320        ONLY:  pmc_childinit, pmc_c_clear_next_array_list,                     &
321               pmc_c_getnextarray, pmc_c_get_2d_index_list, pmc_c_getbuffer,   &
322               pmc_c_putbuffer, pmc_c_setind_and_allocmem,                     &
323               pmc_c_set_dataarray, pmc_set_dataarray_name
324
325    USE pmc_general,                                                           &
326        ONLY:  da_namelen
327
328    USE pmc_handle_communicator,                                               &
329        ONLY:  pmc_get_model_info, pmc_init_model, pmc_is_rootmodel,           &
330               pmc_no_namelist_found, pmc_parent_for_child, m_couplers
331
332    USE pmc_mpi_wrapper,                                                       &
333        ONLY:  pmc_bcast, pmc_recv_from_child, pmc_recv_from_parent,           &
334               pmc_send_to_child, pmc_send_to_parent
335
336    USE pmc_parent,                                                            &
337        ONLY:  pmc_parentinit, pmc_s_clear_next_array_list, pmc_s_fillbuffer,  &
338               pmc_s_getdata_from_buffer, pmc_s_getnextarray,                  &
339               pmc_s_setind_and_allocmem, pmc_s_set_active_data_array,         &
340               pmc_s_set_dataarray, pmc_s_set_2d_index_list
341
342#endif
343
344    USE surface_mod,                                                           &
345        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
346
347    IMPLICIT NONE
348
349#if defined( __parallel )
350#if defined( __mpifh )
351    INCLUDE "mpif.h"
352#endif
353#endif
354
355    PRIVATE
356!
357!-- Constants
358    INTEGER(iwp), PARAMETER ::  child_to_parent = 2   !<
359    INTEGER(iwp), PARAMETER ::  parent_to_child = 1   !<
360!
361!-- Coupler setup
362    INTEGER(iwp), SAVE      ::  comm_world_nesting    !<
363    INTEGER(iwp), SAVE      ::  cpl_id  = 1           !<
364    CHARACTER(LEN=32), SAVE ::  cpl_name              !<
365    INTEGER(iwp), SAVE      ::  cpl_npe_total         !<
366    INTEGER(iwp), SAVE      ::  cpl_parent_id         !<
367!
368!-- Control parameters, will be made input parameters later
369    CHARACTER(LEN=7), SAVE ::  nesting_datatransfer_mode = 'mixed'  !< steering
370                                                         !< parameter for data-
371                                                         !< transfer mode
372    CHARACTER(LEN=8), SAVE ::  nesting_mode = 'two-way'  !< steering parameter
373                                                         !< for 1- or 2-way nesting
374
375    LOGICAL, SAVE ::  nested_run = .FALSE.  !< general switch
376    LOGICAL       ::  rans_mode_parent = .FALSE. !< mode of parent model (.F. - LES mode, .T. - RANS mode)
377
378    REAL(wp), SAVE ::  anterp_relax_length_l = -1.0_wp   !<
379    REAL(wp), SAVE ::  anterp_relax_length_r = -1.0_wp   !<
380    REAL(wp), SAVE ::  anterp_relax_length_s = -1.0_wp   !<
381    REAL(wp), SAVE ::  anterp_relax_length_n = -1.0_wp   !<
382    REAL(wp), SAVE ::  anterp_relax_length_t = -1.0_wp   !<
383!
384!-- Geometry
385    REAL(wp), SAVE                                    ::  area_t             !<
386    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_x            !<
387    REAL(wp), SAVE, DIMENSION(:), ALLOCATABLE, PUBLIC ::  coord_y            !<
388    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_x !<
389    REAL(wp), SAVE, PUBLIC                            ::  lower_left_coord_y !<
390
391!
392!-- Child coarse data arrays
393    INTEGER(iwp), DIMENSION(5),PUBLIC           ::  coarse_bound   !<
394
395    REAL(wp), SAVE                              ::  xexl           !<
396    REAL(wp), SAVE                              ::  xexr           !<
397    REAL(wp), SAVE                              ::  yexs           !<
398    REAL(wp), SAVE                              ::  yexn           !<
399    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_l    !<
400    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_n    !<
401    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_r    !<
402    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_s    !<
403    REAL(wp), SAVE, DIMENSION(:,:), ALLOCATABLE ::  tkefactor_t    !<
404
405    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  dissc !< coarse grid array on child domain - dissipation rate
406    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ec   !<
407    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ptc  !<
408    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uc   !<
409    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vc   !<
410    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc   !<
411    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_c  !<
412    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc  !<
413    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qrc  !<
414    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nrc  !<
415    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !<
416    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc   !<
417    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  nr_partc    !<
418    INTEGER(idp), SAVE, DIMENSION(:,:), ALLOCATABLE, TARGET, PUBLIC ::  part_adrc   !<
419
420    REAL(wp), SAVE, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET ::  chem_spec_c !< coarse grid array on child domain - chemical species
421
422!
423!-- Child interpolation coefficients and child-array indices to be
424!-- precomputed and stored.
425    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ico    !<
426    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  icu    !<
427    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jco    !<
428    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jcv    !<
429    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kco    !<
430    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kcw    !<
431    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xo   !<
432    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xo   !<
433    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1xu   !<
434    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2xu   !<
435    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yo   !<
436    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yo   !<
437    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1yv   !<
438    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2yv   !<
439    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zo   !<
440    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zo   !<
441    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r1zw   !<
442    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:)     ::  r2zw   !<
443!
444!-- Child index arrays and log-ratio arrays for the log-law near-wall
445!-- corrections. These are not truly 3-D arrays but multiple 2-D arrays.
446    INTEGER(iwp), SAVE :: ncorr  !< 4th dimension of the log_ratio-arrays
447    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_l   !<
448    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_n   !<
449    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_r   !<
450    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_u_s   !<
451    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_l   !<
452    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_n   !<
453    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_r   !<
454    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_v_s   !<
455    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_l   !<
456    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_n   !<
457    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_r   !<
458    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::  logc_w_s   !<
459    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_u_l !<
460    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_u_n !<
461    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_u_r !<
462    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_u_s !<
463    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_v_l !<   
464    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_v_n !<
465    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_v_r !<   
466    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_v_s !<
467    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_w_l !<   
468    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_w_n !<
469    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_w_r !<   
470    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::  logc_kbounds_w_s !<       
471    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_l   !<
472    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_n   !<
473    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_r   !<
474    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_u_s   !<
475    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_l   !<
476    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_n   !<
477    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_r   !<
478    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_v_s   !<
479    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_l   !<
480    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_n   !<
481    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_r   !<
482    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:)   ::  logc_ratio_w_s   !<
483!
484!-- Upper bounds for k in anterpolation.
485    INTEGER(iwp), SAVE ::  kctu   !<
486    INTEGER(iwp), SAVE ::  kctw   !<
487!
488!-- Upper bound for k in log-law correction in interpolation.
489    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_l   !<
490    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_n   !<
491    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_r   !<
492    INTEGER(iwp), SAVE ::  nzt_topo_nestbc_s   !<
493!
494!-- Number of ghost nodes in coarse-grid arrays for i and j in anterpolation.
495    INTEGER(iwp), SAVE ::  nhll   !<
496    INTEGER(iwp), SAVE ::  nhlr   !<
497    INTEGER(iwp), SAVE ::  nhls   !<
498    INTEGER(iwp), SAVE ::  nhln   !<
499!
500!-- Spatial under-relaxation coefficients for anterpolation.
501    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  frax   !<
502    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fray   !<
503    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:) ::  fraz   !<
504!
505!-- Child-array indices to be precomputed and stored for anterpolation.
506    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflu   !<
507    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuu   !<
508    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  iflo   !<
509    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  ifuo   !<
510    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflv   !<
511    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuv   !<
512    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jflo   !<
513    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  jfuo   !<
514    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflw   !<
515    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuw   !<
516    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kflo   !<
517    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:) ::  kfuo   !<
518!
519!-- Number of fine-grid nodes inside coarse-grid ij-faces
520!-- to be precomputed for anterpolation.
521    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_u        !<
522    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_v        !<
523    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::  ijfc_s        !<
524    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_w         !<
525    INTEGER(iwp), SAVE, ALLOCATABLE, DIMENSION(:)   ::  kfc_s         !<
526   
527    INTEGER(iwp), DIMENSION(3)          ::  parent_grid_info_int    !<
528    REAL(wp), DIMENSION(7)              ::  parent_grid_info_real   !<
529    REAL(wp), DIMENSION(2)              ::  zmax_coarse             !<
530
531    TYPE coarsegrid_def
532       INTEGER(iwp)                        ::  nx                 !<
533       INTEGER(iwp)                        ::  ny                 !<
534       INTEGER(iwp)                        ::  nz                 !<
535       REAL(wp)                            ::  dx                 !<
536       REAL(wp)                            ::  dy                 !<
537       REAL(wp)                            ::  dz                 !<
538       REAL(wp)                            ::  lower_left_coord_x !<
539       REAL(wp)                            ::  lower_left_coord_y !<
540       REAL(wp)                            ::  xend               !<
541       REAL(wp)                            ::  yend               !<
542       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_x            !<
543       REAL(wp), DIMENSION(:), ALLOCATABLE ::  coord_y            !<
544       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzu                !<
545       REAL(wp), DIMENSION(:), ALLOCATABLE ::  dzw                !<
546       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu                 !<
547       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw                 !<
548    END TYPE coarsegrid_def
549
550    TYPE(coarsegrid_def), SAVE, PUBLIC     ::  cg   !<
551
552!-  Variables for particle coupling
553
554    TYPE, PUBLIC :: childgrid_def
555       INTEGER(iwp)                        ::  nx                   !<
556       INTEGER(iwp)                        ::  ny                   !<
557       INTEGER(iwp)                        ::  nz                   !<
558       REAL(wp)                            ::  dx                   !<
559       REAL(wp)                            ::  dy                   !<
560       REAL(wp)                            ::  dz                   !<
561       REAL(wp)                            ::  lx_coord, lx_coord_b !<
562       REAL(wp)                            ::  rx_coord, rx_coord_b !<
563       REAL(wp)                            ::  sy_coord, sy_coord_b !<
564       REAL(wp)                            ::  ny_coord, ny_coord_b !<
565       REAL(wp)                            ::  uz_coord, uz_coord_b !<
566    END TYPE childgrid_def
567
568    TYPE(childgrid_def), SAVE, ALLOCATABLE, DIMENSION(:), PUBLIC :: childgrid !<
569
570    INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: nr_part  !<
571    INTEGER(idp),ALLOCATABLE,DIMENSION(:,:),PUBLIC,TARGET    :: part_adr !<
572   
573    INTERFACE pmci_boundary_conds
574       MODULE PROCEDURE pmci_boundary_conds
575    END INTERFACE pmci_boundary_conds
576   
577    INTERFACE pmci_check_setting_mismatches
578       MODULE PROCEDURE pmci_check_setting_mismatches
579    END INTERFACE
580
581    INTERFACE pmci_child_initialize
582       MODULE PROCEDURE pmci_child_initialize
583    END INTERFACE
584
585    INTERFACE pmci_synchronize
586       MODULE PROCEDURE pmci_synchronize
587    END INTERFACE
588
589    INTERFACE pmci_datatrans
590       MODULE PROCEDURE pmci_datatrans
591    END INTERFACE pmci_datatrans
592
593    INTERFACE pmci_ensure_nest_mass_conservation
594       MODULE PROCEDURE pmci_ensure_nest_mass_conservation
595    END INTERFACE
596
597    INTERFACE pmci_init
598       MODULE PROCEDURE pmci_init
599    END INTERFACE
600
601    INTERFACE pmci_modelconfiguration
602       MODULE PROCEDURE pmci_modelconfiguration
603    END INTERFACE
604
605    INTERFACE pmci_parent_initialize
606       MODULE PROCEDURE pmci_parent_initialize
607    END INTERFACE
608
609    INTERFACE get_number_of_childs
610       MODULE PROCEDURE get_number_of_childs
611    END  INTERFACE get_number_of_childs
612
613    INTERFACE get_childid
614       MODULE PROCEDURE get_childid
615    END  INTERFACE get_childid
616
617    INTERFACE get_child_edges
618       MODULE PROCEDURE get_child_edges
619    END  INTERFACE get_child_edges
620
621    INTERFACE get_child_gridspacing
622       MODULE PROCEDURE get_child_gridspacing
623    END  INTERFACE get_child_gridspacing
624
625
626    INTERFACE pmci_set_swaplevel
627       MODULE PROCEDURE pmci_set_swaplevel
628    END INTERFACE pmci_set_swaplevel
629
630    PUBLIC anterp_relax_length_l, anterp_relax_length_r,                       &
631           anterp_relax_length_s, anterp_relax_length_n,                       &
632           anterp_relax_length_t, child_to_parent, comm_world_nesting,         &
633           cpl_id, nested_run, nesting_datatransfer_mode, nesting_mode,        &
634           parent_to_child, rans_mode_parent
635
636    PUBLIC pmci_boundary_conds
637    PUBLIC pmci_child_initialize
638    PUBLIC pmci_datatrans
639    PUBLIC pmci_ensure_nest_mass_conservation
640    PUBLIC pmci_init
641    PUBLIC pmci_modelconfiguration
642    PUBLIC pmci_parent_initialize
643    PUBLIC pmci_synchronize
644    PUBLIC pmci_set_swaplevel
645    PUBLIC get_number_of_childs, get_childid, get_child_edges, get_child_gridspacing
646
647
648
649 CONTAINS
650
651
652 SUBROUTINE pmci_init( world_comm )
653
654    USE control_parameters,                                                    &
655        ONLY:  message_string
656
657    IMPLICIT NONE
658
659    INTEGER(iwp), INTENT(OUT) ::  world_comm   !<
660
661#if defined( __parallel )
662
663    INTEGER(iwp)         ::  ierr         !<
664    INTEGER(iwp)         ::  istat        !<
665    INTEGER(iwp)         ::  pmc_status   !<
666
667
668    CALL pmc_init_model( world_comm, nesting_datatransfer_mode, nesting_mode,  &
669                         pmc_status )
670
671    IF ( pmc_status == pmc_no_namelist_found )  THEN
672!
673!--    This is not a nested run
674       world_comm = MPI_COMM_WORLD
675       cpl_id     = 1
676       cpl_name   = ""
677
678       RETURN
679
680    ENDIF
681!
682!-- Check steering parameter values
683    IF ( TRIM( nesting_mode ) /= 'one-way'  .AND.                              &
684         TRIM( nesting_mode ) /= 'two-way'  .AND.                              &
685         TRIM( nesting_mode ) /= 'vertical' )                                  &
686    THEN
687       message_string = 'illegal nesting mode: ' // TRIM( nesting_mode )
688       CALL message( 'pmci_init', 'PA0417', 3, 2, 0, 6, 0 )
689    ENDIF
690
691    IF ( TRIM( nesting_datatransfer_mode ) /= 'cascade'  .AND.                 &
692         TRIM( nesting_datatransfer_mode ) /= 'mixed'    .AND.                 &
693         TRIM( nesting_datatransfer_mode ) /= 'overlap' )                      &
694    THEN
695       message_string = 'illegal nesting datatransfer mode: '                  &
696                        // TRIM( nesting_datatransfer_mode )
697       CALL message( 'pmci_init', 'PA0418', 3, 2, 0, 6, 0 )
698    ENDIF
699!
700!-- Set the general steering switch which tells PALM that its a nested run
701    nested_run = .TRUE.
702!
703!-- Get some variables required by the pmc-interface (and in some cases in the
704!-- PALM code out of the pmci) out of the pmc-core
705    CALL pmc_get_model_info( comm_world_nesting = comm_world_nesting,          &
706                             cpl_id = cpl_id, cpl_parent_id = cpl_parent_id,   &
707                             cpl_name = cpl_name, npe_total = cpl_npe_total,   &
708                             lower_left_x = lower_left_coord_x,                &
709                             lower_left_y = lower_left_coord_y )
710!
711!-- Set the steering switch which tells the models that they are nested (of
712!-- course the root domain (cpl_id = 1) is not nested)
713    IF ( cpl_id >= 2 )  THEN
714       nest_domain = .TRUE.
715       WRITE( coupling_char, '(A2,I2.2)') '_N', cpl_id
716    ENDIF
717
718!
719!-- Message that communicators for nesting are initialized.
720!-- Attention: myid has been set at the end of pmc_init_model in order to
721!-- guarantee that only PE0 of the root domain does the output.
722    CALL location_message( 'finished', .TRUE. )
723!
724!-- Reset myid to its default value
725    myid = 0
726#else
727!
728!-- Nesting cannot be used in serial mode. cpl_id is set to root domain (1)
729!-- because no location messages would be generated otherwise.
730!-- world_comm is given a dummy value to avoid compiler warnings (INTENT(OUT)
731!-- should get an explicit value)
732    cpl_id     = 1
733    nested_run = .FALSE.
734    world_comm = 1
735#endif
736
737 END SUBROUTINE pmci_init
738
739
740
741 SUBROUTINE pmci_modelconfiguration
742
743    IMPLICIT NONE
744
745    INTEGER(iwp) ::  ncpl   !<  number of nest domains
746
747    CALL location_message( 'setup the nested model configuration', .FALSE. )
748    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'start' )
749!
750!-- Compute absolute coordinates for all models
751    CALL pmci_setup_coordinates
752!
753!-- Initialize the child (must be called before pmc_setup_parent)
754    CALL pmci_setup_child
755!
756!-- Initialize PMC parent
757    CALL pmci_setup_parent
758!
759!-- Check for mismatches between settings of master and child variables
760!-- (e.g., all children have to follow the end_time settings of the root master)
761    CALL pmci_check_setting_mismatches
762!
763!-- Set flag file for combine_plot_fields for precessing the nest output data
764    OPEN( 90, FILE='3DNESTING', FORM='FORMATTED' )
765    CALL pmc_get_model_info( ncpl = ncpl )
766    WRITE( 90, '(I2)' )  ncpl
767    CLOSE( 90 )
768
769    CALL cpu_log( log_point_s(79), 'pmci_model_config', 'stop' )
770    CALL location_message( 'finished', .TRUE. )
771
772 END SUBROUTINE pmci_modelconfiguration
773
774
775
776 SUBROUTINE pmci_setup_parent
777
778#if defined( __parallel )
779    IMPLICIT NONE
780
781    CHARACTER(LEN=32) ::  myname
782   
783    INTEGER(iwp) ::  child_id         !<
784    INTEGER(iwp) ::  ierr             !<
785    INTEGER(iwp) ::  i                !<
786    INTEGER(iwp) ::  j                !<
787    INTEGER(iwp) ::  k                !<
788    INTEGER(iwp) ::  m                !<
789    INTEGER(iwp) ::  mid              !<
790    INTEGER(iwp) ::  mm               !<
791    INTEGER(iwp) ::  n = 1            !< running index for chemical species
792    INTEGER(iwp) ::  nest_overlap     !<
793    INTEGER(iwp) ::  nomatch          !<
794    INTEGER(iwp) ::  nx_cl            !<
795    INTEGER(iwp) ::  ny_cl            !<
796    INTEGER(iwp) ::  nz_cl            !<
797
798    INTEGER(iwp), DIMENSION(5) ::  val    !<
799
800
801    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_xl   !<
802    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_xr   !<   
803    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_ys   !<
804    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ch_yn   !<
805    REAL(wp) ::  cl_height        !<
806    REAL(wp) ::  dx_cl            !<
807    REAL(wp) ::  dy_cl            !<
808    REAL(wp) ::  dz_cl            !<
809    REAL(wp) ::  left_limit       !<
810    REAL(wp) ::  north_limit      !<
811    REAL(wp) ::  right_limit      !<
812    REAL(wp) ::  south_limit      !<
813    REAL(wp) ::  xez              !<
814    REAL(wp) ::  yez              !<
815
816    REAL(wp), DIMENSION(5) ::  fval             !<
817
818    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_x   !<
819    REAL(wp), DIMENSION(:), ALLOCATABLE ::  cl_coord_y   !<
820
821!
822!   Initialize the pmc parent
823    CALL pmc_parentinit
824
825!
826!-- Corners of all children of the present parent
827    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN
828       ALLOCATE( ch_xl(1:SIZE( pmc_parent_for_child ) - 1) )
829       ALLOCATE( ch_xr(1:SIZE( pmc_parent_for_child ) - 1) )
830       ALLOCATE( ch_ys(1:SIZE( pmc_parent_for_child ) - 1) )
831       ALLOCATE( ch_yn(1:SIZE( pmc_parent_for_child ) - 1) )
832    ENDIF
833    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) )  THEN
834       ALLOCATE( childgrid(1:SIZE( pmc_parent_for_child ) - 1) )
835    ENDIF
836
837!
838!-- Get coordinates from all children
839    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
840
841       child_id = pmc_parent_for_child(m)
842
843       IF ( myid == 0 )  THEN
844
845          CALL pmc_recv_from_child( child_id, val,  size(val),  0, 123, ierr )
846          CALL pmc_recv_from_child( child_id, fval, size(fval), 0, 124, ierr )
847         
848          nx_cl     = val(1)
849          ny_cl     = val(2)
850          dx_cl     = fval(3)
851          dy_cl     = fval(4)
852          dz_cl     = fval(5)
853          cl_height = fval(1)
854
855          nz_cl = nz
856!
857!--       Find the highest nest level in the coarse grid for the reduced z
858!--       transfer
859          DO  k = 1, nz                 
860             IF ( zw(k) > fval(1) )  THEN
861                nz_cl = k
862                EXIT
863             ENDIF
864          ENDDO
865
866          zmax_coarse = fval(1:2)
867          cl_height   = fval(1)
868
869!   
870!--       Get absolute coordinates from the child
871          ALLOCATE( cl_coord_x(-nbgp:nx_cl+nbgp) )
872          ALLOCATE( cl_coord_y(-nbgp:ny_cl+nbgp) )
873         
874          CALL pmc_recv_from_child( child_id, cl_coord_x, SIZE( cl_coord_x ),  &
875               0, 11, ierr )
876          CALL pmc_recv_from_child( child_id, cl_coord_y, SIZE( cl_coord_y ),  &
877               0, 12, ierr )
878         
879          parent_grid_info_real(1) = lower_left_coord_x
880          parent_grid_info_real(2) = lower_left_coord_y
881          parent_grid_info_real(3) = dx
882          parent_grid_info_real(4) = dy
883          parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx
884          parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy
885          parent_grid_info_real(7) = dz
886
887          parent_grid_info_int(1)  = nx
888          parent_grid_info_int(2)  = ny
889          parent_grid_info_int(3)  = nz_cl
890!
891!--       Check that the child domain matches parent domain.
892          nomatch = 0
893          IF ( nesting_mode == 'vertical' )  THEN
894             right_limit = parent_grid_info_real(5)
895             north_limit = parent_grid_info_real(6)
896             IF ( ( cl_coord_x(nx_cl+1) /= right_limit ) .OR.                  &
897                  ( cl_coord_y(ny_cl+1) /= north_limit ) )  THEN
898                nomatch = 1
899             ENDIF
900          ELSE       
901!
902!--       Check that the child domain is completely inside the parent domain.
903             xez = ( nbgp + 1 ) * dx
904             yez = ( nbgp + 1 ) * dy
905             left_limit  = lower_left_coord_x + xez
906             right_limit = parent_grid_info_real(5) - xez
907             south_limit = lower_left_coord_y + yez
908             north_limit = parent_grid_info_real(6) - yez
909             IF ( ( cl_coord_x(0) < left_limit )        .OR.                   &
910                  ( cl_coord_x(nx_cl+1) > right_limit ) .OR.                   &
911                  ( cl_coord_y(0) < south_limit )       .OR.                   &
912                  ( cl_coord_y(ny_cl+1) > north_limit ) )  THEN
913                nomatch = 1
914             ENDIF
915          ENDIF
916!             
917!--       Child domain must be lower than the parent domain such
918!--       that the top ghost layer of the child grid does not exceed
919!--       the parent domain top boundary.
920
921          IF ( cl_height > zw(nz) ) THEN
922             nomatch = 1
923          ENDIF
924!
925!--       Check that parallel nest domains, if any, do not overlap.
926          nest_overlap = 0
927          IF ( SIZE( pmc_parent_for_child ) - 1 > 0 )  THEN
928             ch_xl(m) = cl_coord_x(-nbgp)
929             ch_xr(m) = cl_coord_x(nx_cl+nbgp)
930             ch_ys(m) = cl_coord_y(-nbgp)
931             ch_yn(m) = cl_coord_y(ny_cl+nbgp)
932
933             IF ( m > 1 )  THEN
934                DO mm = 1, m - 1
935                   mid = pmc_parent_for_child(mm)
936!
937!--                Check Only different nest level
938                   IF (m_couplers(child_id)%parent_id /= m_couplers(mid)%parent_id)  THEN
939                      IF ( ( ch_xl(m) < ch_xr(mm) .OR.                         &
940                             ch_xr(m) > ch_xl(mm) )  .AND.                     &
941                           ( ch_ys(m) < ch_yn(mm) .OR.                         &
942                             ch_yn(m) > ch_ys(mm) ) )  THEN
943                         nest_overlap = 1
944                      ENDIF
945                   ENDIF
946                ENDDO
947             ENDIF
948          ENDIF
949
950          CALL set_child_edge_coords
951
952          DEALLOCATE( cl_coord_x )
953          DEALLOCATE( cl_coord_y )
954
955!
956!--       Send information about operating mode (LES or RANS) to child. This will be
957!--       used to control TKE nesting and setting boundary conditions properly.
958          CALL pmc_send_to_child( child_id, rans_mode, 1, 0, 19, ierr ) 
959!
960!--       Send coarse grid information to child
961          CALL pmc_send_to_child( child_id, parent_grid_info_real,             &
962                                   SIZE( parent_grid_info_real ), 0, 21,       &
963                                   ierr )
964          CALL pmc_send_to_child( child_id, parent_grid_info_int,  3, 0,       &
965                                   22, ierr )
966!
967!--       Send local grid to child
968          CALL pmc_send_to_child( child_id, coord_x, nx+1+2*nbgp, 0, 24,       &
969                                   ierr )
970          CALL pmc_send_to_child( child_id, coord_y, ny+1+2*nbgp, 0, 25,       &
971                                   ierr )
972!
973!--       Also send the dzu-, dzw-, zu- and zw-arrays here
974          CALL pmc_send_to_child( child_id, dzu, nz_cl+1, 0, 26, ierr )
975          CALL pmc_send_to_child( child_id, dzw, nz_cl+1, 0, 27, ierr )
976          CALL pmc_send_to_child( child_id, zu,  nz_cl+2, 0, 28, ierr )
977          CALL pmc_send_to_child( child_id, zw,  nz_cl+2, 0, 29, ierr )
978
979       ENDIF
980
981       CALL MPI_BCAST( nomatch, 1, MPI_INTEGER, 0, comm2d, ierr )
982       IF ( nomatch /= 0 )  THEN
983          WRITE ( message_string, * )  'nested child domain does ',            &
984                                       'not fit into its parent domain'
985          CALL message( 'pmci_setup_parent', 'PA0425', 3, 2, 0, 6, 0 )
986       ENDIF
987 
988       CALL MPI_BCAST( nest_overlap, 1, MPI_INTEGER, 0, comm2d, ierr )
989       IF ( nest_overlap /= 0  .AND.  nesting_mode /= 'vertical' )  THEN
990          WRITE ( message_string, * )  'nested parallel child domains overlap'
991          CALL message( 'pmci_setup_parent', 'PA0426', 3, 2, 0, 6, 0 )
992       ENDIF
993     
994       CALL MPI_BCAST( nz_cl, 1, MPI_INTEGER, 0, comm2d, ierr )
995
996       CALL MPI_BCAST( childgrid(m), STORAGE_SIZE(childgrid(1))/8, MPI_BYTE, 0, comm2d, ierr )
997!
998!--    TO_DO: Klaus: please give a comment what is done here
999       CALL pmci_create_index_list
1000!
1001!--    Include couple arrays into parent content
1002!--    The adresses of the PALM 2D or 3D array (here server coarse grid) which are candidates
1003!--    for coupling are stored once into the pmc context. While data transfer, the array do not
1004!--    have to be specified again
1005
1006       CALL pmc_s_clear_next_array_list
1007       DO  WHILE ( pmc_s_getnextarray( child_id, myname ) )
1008          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
1009             CALL pmci_set_array_pointer( myname, child_id = child_id,         &
1010                                          nz_cl = nz_cl, n = n )
1011             n = n + 1
1012          ELSE
1013             CALL pmci_set_array_pointer( myname, child_id = child_id,         &
1014                                          nz_cl = nz_cl )
1015          ENDIF
1016       ENDDO
1017
1018       CALL pmc_s_setind_and_allocmem( child_id )
1019    ENDDO
1020
1021    IF ( ( SIZE( pmc_parent_for_child ) - 1 > 0 ) .AND. myid == 0 )  THEN
1022       DEALLOCATE( ch_xl )
1023       DEALLOCATE( ch_xr )
1024       DEALLOCATE( ch_ys )
1025       DEALLOCATE( ch_yn )
1026    ENDIF
1027
1028 CONTAINS
1029
1030
1031   SUBROUTINE pmci_create_index_list
1032
1033       IMPLICIT NONE
1034
1035       INTEGER(iwp) ::  i                  !<
1036       INTEGER(iwp) ::  ic                 !<
1037       INTEGER(iwp) ::  ierr               !<
1038       INTEGER(iwp) ::  j                  !<
1039       INTEGER(iwp) ::  k                  !<
1040       INTEGER(iwp) ::  m                  !<
1041       INTEGER(iwp) ::  n                  !<
1042       INTEGER(iwp) ::  npx                !<
1043       INTEGER(iwp) ::  npy                !<
1044       INTEGER(iwp) ::  nrx                !<
1045       INTEGER(iwp) ::  nry                !<
1046       INTEGER(iwp) ::  px                 !<
1047       INTEGER(iwp) ::  py                 !<
1048       INTEGER(iwp) ::  parent_pe          !<
1049
1050       INTEGER(iwp), DIMENSION(2) ::  scoord             !<
1051       INTEGER(iwp), DIMENSION(2) ::  size_of_array      !<
1052
1053       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  coarse_bound_all   !<
1054       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  ::  index_list         !<
1055
1056       IF ( myid == 0 )  THEN
1057!         
1058!--       TO_DO: Klaus: give more specific comment what size_of_array stands for
1059          CALL pmc_recv_from_child( child_id, size_of_array, 2, 0, 40, ierr )
1060          ALLOCATE( coarse_bound_all(size_of_array(1),size_of_array(2)) )
1061          CALL pmc_recv_from_child( child_id, coarse_bound_all,                &
1062                                    SIZE( coarse_bound_all ), 0, 41, ierr )
1063!
1064!--       Compute size of index_list.
1065          ic = 0
1066          DO  k = 1, size_of_array(2)
1067             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
1068                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
1069                   ic = ic + 1
1070                ENDDO
1071             ENDDO
1072          ENDDO
1073
1074          ALLOCATE( index_list(6,ic) )
1075
1076          CALL MPI_COMM_SIZE( comm1dx, npx, ierr )
1077          CALL MPI_COMM_SIZE( comm1dy, npy, ierr )
1078!
1079!--       The +1 in index is because PALM starts with nx=0
1080          nrx = nxr - nxl + 1
1081          nry = nyn - nys + 1
1082          ic  = 0
1083!
1084!--       Loop over all children PEs
1085          DO  k = 1, size_of_array(2)
1086!
1087!--          Area along y required by actual child PE
1088             DO  j = coarse_bound_all(3,k), coarse_bound_all(4,k)
1089!
1090!--             Area along x required by actual child PE
1091                DO  i = coarse_bound_all(1,k), coarse_bound_all(2,k)
1092
1093                   px = i / nrx
1094                   py = j / nry
1095                   scoord(1) = px
1096                   scoord(2) = py
1097                   CALL MPI_CART_RANK( comm2d, scoord, parent_pe, ierr )
1098                 
1099                   ic = ic + 1
1100!
1101!--                First index in parent array
1102                   index_list(1,ic) = i - ( px * nrx ) + 1 + nbgp
1103!
1104!--                Second index in parent array
1105                   index_list(2,ic) = j - ( py * nry ) + 1 + nbgp
1106!
1107!--                x index of child coarse grid
1108                   index_list(3,ic) = i - coarse_bound_all(1,k) + 1
1109!
1110!--                y index of child coarse grid
1111                   index_list(4,ic) = j - coarse_bound_all(3,k) + 1
1112!
1113!--                PE number of child
1114                   index_list(5,ic) = k - 1
1115!
1116!--                PE number of parent
1117                   index_list(6,ic) = parent_pe
1118
1119                ENDDO
1120             ENDDO
1121          ENDDO
1122!
1123!--       TO_DO: Klaus: comment what is done here
1124          CALL pmc_s_set_2d_index_list( child_id, index_list(:,1:ic) )
1125
1126       ELSE
1127!
1128!--       TO_DO: Klaus: comment why this dummy allocation is required
1129          ALLOCATE( index_list(6,1) )
1130          CALL pmc_s_set_2d_index_list( child_id, index_list )
1131       ENDIF
1132
1133       DEALLOCATE(index_list)
1134
1135     END SUBROUTINE pmci_create_index_list
1136
1137     SUBROUTINE set_child_edge_coords
1138        IMPLICIT  NONE
1139
1140        INTEGER(iwp) :: nbgp_lpm = 1
1141
1142        nbgp_lpm = min(nbgp_lpm, nbgp)
1143
1144        childgrid(m)%nx = nx_cl
1145        childgrid(m)%ny = ny_cl
1146        childgrid(m)%nz = nz_cl
1147        childgrid(m)%dx = dx_cl
1148        childgrid(m)%dy = dy_cl
1149        childgrid(m)%dz = dz_cl
1150
1151        childgrid(m)%lx_coord   = cl_coord_x(0)
1152        childgrid(m)%lx_coord_b = cl_coord_x(-nbgp_lpm)
1153        childgrid(m)%rx_coord   = cl_coord_x(nx_cl)+dx_cl
1154        childgrid(m)%rx_coord_b = cl_coord_x(nx_cl+nbgp_lpm)+dx_cl
1155        childgrid(m)%sy_coord   = cl_coord_y(0)
1156        childgrid(m)%sy_coord_b = cl_coord_y(-nbgp_lpm)
1157        childgrid(m)%ny_coord   = cl_coord_y(ny_cl)+dy_cl
1158        childgrid(m)%ny_coord_b = cl_coord_y(ny_cl+nbgp_lpm)+dy_cl
1159        childgrid(m)%uz_coord   = zmax_coarse(2)
1160        childgrid(m)%uz_coord_b = zmax_coarse(1)
1161
1162!         WRITE(9,*)                 'edge coordinates for child id ',child_id,m
1163!         WRITE(9,*)                 'Number of Boundray cells lpm  ',nbgp_lpm
1164!         WRITE(9,'(a,3i7,2f10.2)') ' model size                    ', nx_cl, ny_cl, nz_cl, dx_cl, dy_cl
1165!          WRITE(9,'(a,5f10.2)')     ' model edge                    ', childgrid(m)%lx_coord,  &
1166!                                childgrid(m)%rx_coord, childgrid(m)%sy_coord,                  &
1167!                                childgrid(m)%ny_coord,childgrid(m)%uz_coord
1168!          WRITE(9,'(a,4f10.2)')     ' model edge with Boundary      ', childgrid(m)%lx_coord_b,&
1169!                                childgrid(m)%rx_coord_b, childgrid(m)%sy_coord_b,              &
1170!                                childgrid(m)%ny_coord_b
1171
1172     END SUBROUTINE set_child_edge_coords
1173
1174#endif
1175 END SUBROUTINE pmci_setup_parent
1176
1177
1178
1179 SUBROUTINE pmci_setup_child
1180
1181
1182#if defined( __parallel )
1183    IMPLICIT NONE
1184
1185    CHARACTER(LEN=da_namelen) ::  myname     !<
1186
1187    INTEGER(iwp) ::  i          !<
1188    INTEGER(iwp) ::  ierr       !<
1189    INTEGER(iwp) ::  icl        !<
1190    INTEGER(iwp) ::  icr        !<
1191    INTEGER(iwp) ::  j          !<
1192    INTEGER(iwp) ::  jcn        !<
1193    INTEGER(iwp) ::  jcs        !<
1194    INTEGER(iwp) ::  n          !< running index for number of chemical species
1195
1196    INTEGER(iwp), DIMENSION(5) ::  val        !<
1197   
1198    REAL(wp) ::  xcs        !<
1199    REAL(wp) ::  xce        !<
1200    REAL(wp) ::  ycs        !<
1201    REAL(wp) ::  yce        !<
1202
1203    REAL(wp), DIMENSION(5) ::  fval       !<
1204                                             
1205!
1206!-- Child setup
1207!-- Root model does not have a parent and is not a child, therefore no child setup on root model
1208
1209    IF ( .NOT. pmc_is_rootmodel() )  THEN
1210
1211       CALL pmc_childinit
1212!
1213!--    Here AND ONLY HERE the arrays are defined, which actualy will be
1214!--    exchanged between child and parent.
1215!--    If a variable is removed, it only has to be removed from here.
1216!--    Please check, if the arrays are in the list of POSSIBLE exchange arrays
1217!--    in subroutines:
1218!--    pmci_set_array_pointer (for parent arrays)
1219!--    pmci_create_child_arrays (for child arrays)
1220       CALL pmc_set_dataarray_name( 'coarse', 'u'  ,'fine', 'u',  ierr )
1221       CALL pmc_set_dataarray_name( 'coarse', 'v'  ,'fine', 'v',  ierr )
1222       CALL pmc_set_dataarray_name( 'coarse', 'w'  ,'fine', 'w',  ierr )
1223!
1224!--    Set data array name for TKE. Please note, nesting of TKE is actually
1225!--    only done if both parent and child are in LES or in RANS mode. Due to
1226!--    design of model coupler, however, data array names must be already
1227!--    available at this point.
1228       CALL pmc_set_dataarray_name( 'coarse', 'e'  ,'fine', 'e',  ierr )
1229!
1230!--    Nesting of dissipation rate only if both parent and child are in RANS
1231!--    mode and TKE-epsilo closure is applied. Please so also comment for TKE
1232!--    above.
1233       CALL pmc_set_dataarray_name( 'coarse', 'diss'  ,'fine', 'diss',  ierr )
1234
1235       IF ( .NOT. neutral )  THEN
1236          CALL pmc_set_dataarray_name( 'coarse', 'pt' ,'fine', 'pt', ierr )
1237       ENDIF
1238
1239       IF ( humidity )  THEN
1240
1241          CALL pmc_set_dataarray_name( 'coarse', 'q'  ,'fine', 'q',  ierr )
1242
1243          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
1244            CALL pmc_set_dataarray_name( 'coarse', 'qc'  ,'fine', 'qc',  ierr ) 
1245            CALL pmc_set_dataarray_name( 'coarse', 'nc'  ,'fine', 'nc',  ierr ) 
1246          ENDIF
1247
1248          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
1249             CALL pmc_set_dataarray_name( 'coarse', 'qr'  ,'fine', 'qr',  ierr )
1250             CALL pmc_set_dataarray_name( 'coarse', 'nr'  ,'fine', 'nr',  ierr ) 
1251          ENDIF
1252     
1253       ENDIF
1254
1255       IF ( passive_scalar )  THEN
1256          CALL pmc_set_dataarray_name( 'coarse', 's'  ,'fine', 's',  ierr )
1257       ENDIF
1258
1259       IF( particle_advection )  THEN
1260          CALL pmc_set_dataarray_name( 'coarse', 'nr_part'  ,'fine',           &
1261               'nr_part',  ierr )
1262          CALL pmc_set_dataarray_name( 'coarse', 'part_adr'  ,'fine',          &
1263               'part_adr',  ierr )
1264       ENDIF
1265       
1266       IF ( air_chemistry )  THEN
1267          DO  n = 1, nspec
1268             CALL pmc_set_dataarray_name( 'coarse',                            &
1269                                          'chem_' //                           &
1270                                          TRIM( chem_species(n)%name ),        &
1271                                         'fine',                               &
1272                                          'chem_' //                           &
1273                                          TRIM( chem_species(n)%name ),        &
1274                                          ierr )
1275          ENDDO
1276       ENDIF
1277
1278       CALL pmc_set_dataarray_name( lastentry = .TRUE. )
1279!
1280!--    Send grid to parent
1281       val(1)  = nx
1282       val(2)  = ny
1283       val(3)  = nz
1284       val(4)  = dx
1285       val(5)  = dy
1286       fval(1) = zw(nzt+1)
1287       fval(2) = zw(nzt)
1288       fval(3) = dx
1289       fval(4) = dy
1290       fval(5) = dz
1291
1292       IF ( myid == 0 )  THEN
1293
1294          CALL pmc_send_to_parent( val, SIZE( val ), 0, 123, ierr )
1295          CALL pmc_send_to_parent( fval, SIZE( fval ), 0, 124, ierr )
1296          CALL pmc_send_to_parent( coord_x, nx + 1 + 2 * nbgp, 0, 11, ierr )
1297          CALL pmc_send_to_parent( coord_y, ny + 1 + 2 * nbgp, 0, 12, ierr )
1298
1299          CALL pmc_recv_from_parent( rans_mode_parent, 1, 0, 19, ierr )
1300!
1301!
1302!--       Receive Coarse grid information.
1303          CALL pmc_recv_from_parent( parent_grid_info_real,                    &
1304                                     SIZE(parent_grid_info_real), 0, 21, ierr )
1305          CALL pmc_recv_from_parent( parent_grid_info_int,  3, 0, 22, ierr )
1306!
1307!--        Debug-printouts - keep them
1308!           WRITE(0,*) 'Coarse grid from parent '
1309!           WRITE(0,*) 'startx_tot    = ',parent_grid_info_real(1)
1310!           WRITE(0,*) 'starty_tot    = ',parent_grid_info_real(2)
1311!           WRITE(0,*) 'endx_tot      = ',parent_grid_info_real(5)
1312!           WRITE(0,*) 'endy_tot      = ',parent_grid_info_real(6)
1313!           WRITE(0,*) 'dx            = ',parent_grid_info_real(3)
1314!           WRITE(0,*) 'dy            = ',parent_grid_info_real(4)
1315!           WRITE(0,*) 'dz            = ',parent_grid_info_real(7)
1316!           WRITE(0,*) 'nx_coarse     = ',parent_grid_info_int(1)
1317!           WRITE(0,*) 'ny_coarse     = ',parent_grid_info_int(2)
1318!           WRITE(0,*) 'nz_coarse     = ',parent_grid_info_int(3)
1319       ENDIF
1320
1321       CALL MPI_BCAST( parent_grid_info_real, SIZE(parent_grid_info_real),     &
1322                       MPI_REAL, 0, comm2d, ierr )
1323       CALL MPI_BCAST( parent_grid_info_int, 3, MPI_INTEGER, 0, comm2d, ierr )
1324
1325       cg%dx = parent_grid_info_real(3)
1326       cg%dy = parent_grid_info_real(4)
1327       cg%dz = parent_grid_info_real(7)
1328       cg%nx = parent_grid_info_int(1)
1329       cg%ny = parent_grid_info_int(2)
1330       cg%nz = parent_grid_info_int(3)
1331!
1332!--    Get parent coordinates on coarse grid
1333       ALLOCATE( cg%coord_x(-nbgp:cg%nx+nbgp) )
1334       ALLOCATE( cg%coord_y(-nbgp:cg%ny+nbgp) )
1335     
1336       ALLOCATE( cg%dzu(1:cg%nz+1) )
1337       ALLOCATE( cg%dzw(1:cg%nz+1) )
1338       ALLOCATE( cg%zu(0:cg%nz+1) )
1339       ALLOCATE( cg%zw(0:cg%nz+1) )
1340!
1341!--    Get coarse grid coordinates and values of the z-direction from the parent
1342       IF ( myid == 0)  THEN
1343
1344          CALL pmc_recv_from_parent( cg%coord_x, cg%nx+1+2*nbgp, 0, 24, ierr )
1345          CALL pmc_recv_from_parent( cg%coord_y, cg%ny+1+2*nbgp, 0, 25, ierr )
1346          CALL pmc_recv_from_parent( cg%dzu, cg%nz + 1, 0, 26, ierr )
1347          CALL pmc_recv_from_parent( cg%dzw, cg%nz + 1, 0, 27, ierr )
1348          CALL pmc_recv_from_parent( cg%zu, cg%nz + 2, 0, 28, ierr )
1349          CALL pmc_recv_from_parent( cg%zw, cg%nz + 2, 0, 29, ierr )
1350
1351       ENDIF
1352!
1353!--    Broadcast this information
1354       CALL MPI_BCAST( cg%coord_x, cg%nx+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1355       CALL MPI_BCAST( cg%coord_y, cg%ny+1+2*nbgp, MPI_REAL, 0, comm2d, ierr )
1356       CALL MPI_BCAST( cg%dzu, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
1357       CALL MPI_BCAST( cg%dzw, cg%nz+1, MPI_REAL, 0, comm2d, ierr )
1358       CALL MPI_BCAST( cg%zu, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
1359       CALL MPI_BCAST( cg%zw, cg%nz+2,  MPI_REAL, 0, comm2d, ierr )
1360       CALL MPI_BCAST( rans_mode_parent, 1, MPI_LOGICAL, 0, comm2d, ierr )
1361 
1362!
1363!--    Find the index bounds for the nest domain in the coarse-grid index space
1364       CALL pmci_map_fine_to_coarse_grid
1365!
1366!--    TO_DO: Klaus give a comment what is happening here
1367       CALL pmc_c_get_2d_index_list
1368!
1369!--    Include couple arrays into child content
1370!--    TO_DO: Klaus: better explain the above comment (what is child content?)
1371       CALL  pmc_c_clear_next_array_list
1372
1373       n = 1
1374       DO  WHILE ( pmc_c_getnextarray( myname ) )
1375!--       Note that cg%nz is not the original nz of parent, but the highest
1376!--       parent-grid level needed for nesting.
1377!--       Please note, in case of chemical species an additional parameter
1378!--       need to be passed, which is required to set the pointer correctly
1379!--       to the chemical-species data structure. Hence, first check if current
1380!--       variable is a chemical species. If so, pass index id of respective
1381!--       species and increment this subsequently.
1382          IF ( INDEX( TRIM( myname ), 'chem_' ) /= 0 )  THEN             
1383             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz, n )
1384             n = n + 1
1385          ELSE
1386             CALL pmci_create_child_arrays ( myname, icl, icr, jcs, jcn, cg%nz )
1387          ENDIF
1388       ENDDO
1389       CALL pmc_c_setind_and_allocmem
1390!
1391!--    Precompute interpolation coefficients and child-array indices
1392       CALL pmci_init_interp_tril
1393!
1394!--    Precompute the log-law correction index- and ratio-arrays
1395       IF ( constant_flux_layer )  THEN
1396          CALL pmci_init_loglaw_correction
1397       ENDIF
1398!
1399!--    Define the SGS-TKE scaling factor based on the grid-spacing ratio. Only
1400!--    if both parent and child are in LES mode or in RANS mode.
1401!--    Please note, in case parent and child are in RANS mode, TKE weighting
1402!--    factor is simply one.
1403       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
1404            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
1405               .NOT. constant_diffusion ) )  CALL pmci_init_tkefactor
1406!
1407!--    Two-way coupling for general and vertical nesting.
1408!--    Precompute the index arrays and relaxation functions for the
1409!--    anterpolation
1410       IF ( TRIM( nesting_mode ) == 'two-way' .OR.                             &
1411                  nesting_mode == 'vertical' )  THEN
1412          CALL pmci_init_anterp_tophat
1413       ENDIF
1414!
1415!--    Finally, compute the total area of the top-boundary face of the domain.
1416!--    This is needed in the pmc_ensure_nest_mass_conservation     
1417       area_t = ( nx + 1 ) * (ny + 1 ) * dx * dy
1418
1419    ENDIF
1420
1421 CONTAINS
1422
1423
1424    SUBROUTINE pmci_map_fine_to_coarse_grid
1425!
1426!--    Determine index bounds of interpolation/anterpolation area in the coarse
1427!--    grid index space
1428       IMPLICIT NONE
1429
1430       INTEGER(iwp), DIMENSION(5,numprocs) ::  coarse_bound_all   !<
1431       INTEGER(iwp), DIMENSION(2)          ::  size_of_array      !<
1432                                             
1433       REAL(wp) ::  loffset     !<
1434       REAL(wp) ::  noffset     !<
1435       REAL(wp) ::  roffset     !<
1436       REAL(wp) ::  soffset     !<
1437
1438!
1439!--    If the fine- and coarse grid nodes do not match:
1440       loffset = MOD( coord_x(nxl), cg%dx )
1441       xexl    = cg%dx + loffset
1442!
1443!--    This is needed in the anterpolation phase
1444       nhll = CEILING( xexl / cg%dx )
1445       xcs  = coord_x(nxl) - xexl
1446       DO  i = 0, cg%nx
1447          IF ( cg%coord_x(i) > xcs )  THEN
1448             icl = MAX( -1, i-1 )
1449             EXIT
1450          ENDIF
1451       ENDDO
1452!
1453!--    If the fine- and coarse grid nodes do not match
1454       roffset = MOD( coord_x(nxr+1), cg%dx )
1455       xexr    = cg%dx + roffset
1456!
1457!--    This is needed in the anterpolation phase
1458       nhlr = CEILING( xexr / cg%dx )
1459       xce  = coord_x(nxr+1) + xexr
1460!--    One "extra" layer is taken behind the right boundary
1461!--    because it may be needed in cases of non-integer grid-spacing ratio
1462       DO  i = cg%nx, 0 , -1
1463          IF ( cg%coord_x(i) < xce )  THEN
1464             icr = MIN( cg%nx+1, i+1 )
1465             EXIT
1466          ENDIF
1467       ENDDO
1468!
1469!--    If the fine- and coarse grid nodes do not match
1470       soffset = MOD( coord_y(nys), cg%dy )
1471       yexs    = cg%dy + soffset
1472!
1473!--    This is needed in the anterpolation phase
1474       nhls = CEILING( yexs / cg%dy )
1475       ycs  = coord_y(nys) - yexs
1476       DO  j = 0, cg%ny
1477          IF ( cg%coord_y(j) > ycs )  THEN
1478             jcs = MAX( -nbgp, j-1 )
1479             EXIT
1480          ENDIF
1481       ENDDO
1482!
1483!--    If the fine- and coarse grid nodes do not match
1484       noffset = MOD( coord_y(nyn+1), cg%dy )
1485       yexn    = cg%dy + noffset
1486!
1487!--    This is needed in the anterpolation phase
1488       nhln = CEILING( yexn / cg%dy )
1489       yce  = coord_y(nyn+1) + yexn
1490!--    One "extra" layer is taken behind the north boundary
1491!--    because it may be needed in cases of non-integer grid-spacing ratio
1492       DO  j = cg%ny, 0, -1
1493          IF ( cg%coord_y(j) < yce )  THEN
1494             jcn = MIN( cg%ny + nbgp, j+1 )
1495             EXIT
1496          ENDIF
1497       ENDDO
1498
1499       coarse_bound(1) = icl
1500       coarse_bound(2) = icr
1501       coarse_bound(3) = jcs
1502       coarse_bound(4) = jcn
1503       coarse_bound(5) = myid
1504!
1505!--    Note that MPI_Gather receives data from all processes in the rank order
1506!--    TO_DO: refer to the line where this fact becomes important
1507       CALL MPI_GATHER( coarse_bound, 5, MPI_INTEGER, coarse_bound_all, 5,     &
1508                        MPI_INTEGER, 0, comm2d, ierr )
1509
1510       IF ( myid == 0 )  THEN
1511          size_of_array(1) = SIZE( coarse_bound_all, 1 )
1512          size_of_array(2) = SIZE( coarse_bound_all, 2 )
1513          CALL pmc_send_to_parent( size_of_array, 2, 0, 40, ierr )
1514          CALL pmc_send_to_parent( coarse_bound_all, SIZE( coarse_bound_all ), &
1515                                   0, 41, ierr )
1516       ENDIF
1517
1518    END SUBROUTINE pmci_map_fine_to_coarse_grid
1519
1520
1521
1522    SUBROUTINE pmci_init_interp_tril
1523!
1524!--    Precomputation of the interpolation coefficients and child-array indices
1525!--    to be used by the interpolation routines interp_tril_lr, interp_tril_ns
1526!--    and interp_tril_t.
1527
1528       IMPLICIT NONE
1529
1530       INTEGER(iwp) ::  i       !<
1531       INTEGER(iwp) ::  i1      !<
1532       INTEGER(iwp) ::  j       !<
1533       INTEGER(iwp) ::  j1      !<
1534       INTEGER(iwp) ::  k       !<
1535       INTEGER(iwp) ::  kc      !<
1536       INTEGER(iwp) ::  kdzo    !<
1537       INTEGER(iwp) ::  kdzw    !<       
1538
1539       REAL(wp) ::  xb          !<
1540       REAL(wp) ::  xcsu        !<
1541       REAL(wp) ::  xfso        !<
1542       REAL(wp) ::  xcso        !<
1543       REAL(wp) ::  xfsu        !<
1544       REAL(wp) ::  yb          !<
1545       REAL(wp) ::  ycso        !<
1546       REAL(wp) ::  ycsv        !<
1547       REAL(wp) ::  yfso        !<
1548       REAL(wp) ::  yfsv        !<
1549       REAL(wp) ::  zcso        !<
1550       REAL(wp) ::  zcsw        !<
1551       REAL(wp) ::  zfso        !<
1552       REAL(wp) ::  zfsw        !<
1553     
1554
1555       xb = nxl * dx
1556       yb = nys * dy
1557     
1558       ALLOCATE( icu(nxlg:nxrg) )
1559       ALLOCATE( ico(nxlg:nxrg) )
1560       ALLOCATE( jcv(nysg:nyng) )
1561       ALLOCATE( jco(nysg:nyng) )
1562       ALLOCATE( kcw(nzb:nzt+1) )
1563       ALLOCATE( kco(nzb:nzt+1) )
1564       ALLOCATE( r1xu(nxlg:nxrg) )
1565       ALLOCATE( r2xu(nxlg:nxrg) )
1566       ALLOCATE( r1xo(nxlg:nxrg) )
1567       ALLOCATE( r2xo(nxlg:nxrg) )
1568       ALLOCATE( r1yv(nysg:nyng) )
1569       ALLOCATE( r2yv(nysg:nyng) )
1570       ALLOCATE( r1yo(nysg:nyng) )
1571       ALLOCATE( r2yo(nysg:nyng) )
1572       ALLOCATE( r1zw(nzb:nzt+1) )
1573       ALLOCATE( r2zw(nzb:nzt+1) )
1574       ALLOCATE( r1zo(nzb:nzt+1) )
1575       ALLOCATE( r2zo(nzb:nzt+1) )
1576!
1577!--    Note that the node coordinates xfs... and xcs... are relative to the
1578!--    lower-left-bottom corner of the fc-array, not the actual child domain
1579!--    corner
1580       DO  i = nxlg, nxrg
1581          xfsu    = coord_x(i) - ( lower_left_coord_x + xb - xexl )
1582          xfso    = coord_x(i) + 0.5_wp * dx - ( lower_left_coord_x + xb - xexl )
1583          icu(i)  = icl + FLOOR( xfsu / cg%dx )
1584          ico(i)  = icl + FLOOR( ( xfso - 0.5_wp * cg%dx ) / cg%dx )
1585          xcsu    = ( icu(i) - icl ) * cg%dx
1586          xcso    = ( ico(i) - icl ) * cg%dx + 0.5_wp * cg%dx
1587          r2xu(i) = ( xfsu - xcsu ) / cg%dx
1588          r2xo(i) = ( xfso - xcso ) / cg%dx
1589          r1xu(i) = 1.0_wp - r2xu(i)
1590          r1xo(i) = 1.0_wp - r2xo(i)
1591       ENDDO
1592
1593       DO  j = nysg, nyng
1594          yfsv    = coord_y(j) - ( lower_left_coord_y + yb - yexs )
1595          yfso    = coord_y(j) + 0.5_wp * dy - ( lower_left_coord_y + yb - yexs )
1596          jcv(j)  = jcs + FLOOR( yfsv / cg%dy )
1597          jco(j)  = jcs + FLOOR( ( yfso -0.5_wp * cg%dy ) / cg%dy )
1598          ycsv    = ( jcv(j) - jcs ) * cg%dy
1599          ycso    = ( jco(j) - jcs ) * cg%dy + 0.5_wp * cg%dy
1600          r2yv(j) = ( yfsv - ycsv ) / cg%dy
1601          r2yo(j) = ( yfso - ycso ) / cg%dy
1602          r1yv(j) = 1.0_wp - r2yv(j)
1603          r1yo(j) = 1.0_wp - r2yo(j)
1604       ENDDO
1605
1606       DO  k = nzb, nzt + 1
1607          zfsw = zw(k)
1608          zfso = zu(k)
1609
1610          DO kc = 0, cg%nz+1
1611             IF ( cg%zw(kc) > zfsw )  EXIT
1612          ENDDO
1613          kcw(k) = kc - 1
1614         
1615          DO kc = 0, cg%nz+1
1616             IF ( cg%zu(kc) > zfso )  EXIT
1617          ENDDO
1618          kco(k) = kc - 1
1619
1620          zcsw    = cg%zw(kcw(k))
1621          zcso    = cg%zu(kco(k))
1622          kdzw    = MIN( kcw(k)+1, cg%nz+1 )
1623          kdzo    = MIN( kco(k)+1, cg%nz+1 )
1624          r2zw(k) = ( zfsw - zcsw ) / cg%dzw(kdzw)
1625          r2zo(k) = ( zfso - zcso ) / cg%dzu(kdzo)
1626          r1zw(k) = 1.0_wp - r2zw(k)
1627          r1zo(k) = 1.0_wp - r2zo(k)
1628       ENDDO
1629     
1630    END SUBROUTINE pmci_init_interp_tril
1631
1632
1633
1634    SUBROUTINE pmci_init_loglaw_correction
1635!
1636!--    Precomputation of the index and log-ratio arrays for the log-law
1637!--    corrections for near-wall nodes after the nest-BC interpolation.
1638!--    These are used by the interpolation routines interp_tril_lr and
1639!--    interp_tril_ns.
1640
1641       IMPLICIT NONE
1642
1643       INTEGER(iwp) ::  direction      !< Wall normal index: 1=k, 2=j, 3=i.
1644       INTEGER(iwp) ::  dum            !< dummy value for reduce operation
1645       INTEGER(iwp) ::  i              !<
1646       INTEGER(iwp) ::  icorr          !<
1647       INTEGER(iwp) ::  ierr           !< MPI status
1648       INTEGER(iwp) ::  inc            !< Wall outward-normal index increment -1
1649                                       !< or 1, for direction=1, inc=1 always
1650       INTEGER(iwp) ::  iw             !<
1651       INTEGER(iwp) ::  j              !<
1652       INTEGER(iwp) ::  jcorr          !<
1653       INTEGER(iwp) ::  jw             !<
1654       INTEGER(iwp) ::  k              !<
1655       INTEGER(iwp) ::  k_wall_u_ji    !< topography top index on u-grid
1656       INTEGER(iwp) ::  k_wall_u_ji_p  !< topography top index on u-grid
1657       INTEGER(iwp) ::  k_wall_u_ji_m  !< topography top index on u-grid
1658       INTEGER(iwp) ::  k_wall_v_ji    !< topography top index on v-grid
1659       INTEGER(iwp) ::  k_wall_v_ji_p  !< topography top index on v-grid
1660       INTEGER(iwp) ::  k_wall_v_ji_m  !< topography top index on v-grid
1661       INTEGER(iwp) ::  k_wall_w_ji    !< topography top index on w-grid
1662       INTEGER(iwp) ::  k_wall_w_ji_p  !< topography top index on w-grid
1663       INTEGER(iwp) ::  k_wall_w_ji_m  !< topography top index on w-grid
1664       INTEGER(iwp) ::  kb             !<
1665       INTEGER(iwp) ::  kcorr          !<
1666       INTEGER(iwp) ::  lc             !<
1667       INTEGER(iwp) ::  m              !< Running index for surface data type
1668       INTEGER(iwp) ::  ni             !<
1669       INTEGER(iwp) ::  nj             !<
1670       INTEGER(iwp) ::  nk             !<
1671       INTEGER(iwp) ::  nzt_topo_max   !<
1672       INTEGER(iwp) ::  wall_index     !<  Index of the wall-node coordinate
1673
1674       REAL(wp)     ::  z0_topo      !<  roughness at vertical walls
1675       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  lcr   !<
1676
1677!
1678!--    First determine the maximum k-index needed for the near-wall corrections.
1679!--    This maximum is individual for each boundary to minimize the storage
1680!--    requirements and to minimize the corresponding loop k-range in the
1681!--    interpolation routines.
1682       nzt_topo_nestbc_l = nzb
1683       IF ( nest_bound_l )  THEN
1684          DO  i = nxl-1, nxl
1685             DO  j = nys, nyn
1686!
1687!--             Concept need to be reconsidered for 3D-topography
1688!--             Determine largest topography index on scalar grid
1689                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
1690                                    get_topography_top_index_ji( j, i, 's' ) )
1691!
1692!--             Determine largest topography index on u grid
1693                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
1694                                    get_topography_top_index_ji( j, i, 'u' ) )
1695!
1696!--             Determine largest topography index on v grid
1697                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
1698                                    get_topography_top_index_ji( j, i, 'v' ) )
1699!
1700!--             Determine largest topography index on w grid
1701                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
1702                                    get_topography_top_index_ji( j, i, 'w' ) )
1703             ENDDO
1704          ENDDO
1705          nzt_topo_nestbc_l = nzt_topo_nestbc_l + 1
1706       ENDIF
1707     
1708       nzt_topo_nestbc_r = nzb
1709       IF ( nest_bound_r )  THEN
1710          i = nxr + 1
1711          DO  j = nys, nyn
1712!
1713!--             Concept need to be reconsidered for 3D-topography
1714!--             Determine largest topography index on scalar grid
1715                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
1716                                    get_topography_top_index_ji( j, i, 's' ) )
1717!
1718!--             Determine largest topography index on u grid
1719                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
1720                                    get_topography_top_index_ji( j, i, 'u' ) )
1721!
1722!--             Determine largest topography index on v grid
1723                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
1724                                    get_topography_top_index_ji( j, i, 'v' ) )
1725!
1726!--             Determine largest topography index on w grid
1727                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
1728                                    get_topography_top_index_ji( j, i, 'w' ) )
1729          ENDDO
1730          nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
1731       ENDIF
1732
1733       nzt_topo_nestbc_s = nzb
1734       IF ( nest_bound_s )  THEN
1735          DO  j = nys-1, nys
1736             DO  i = nxl, nxr
1737!
1738!--             Concept need to be reconsidered for 3D-topography
1739!--             Determine largest topography index on scalar grid
1740                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
1741                                    get_topography_top_index_ji( j, i, 's' ) )
1742!
1743!--             Determine largest topography index on u grid
1744                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
1745                                    get_topography_top_index_ji( j, i, 'u' ) )
1746!
1747!--             Determine largest topography index on v grid
1748                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
1749                                    get_topography_top_index_ji( j, i, 'v' ) )
1750!
1751!--             Determine largest topography index on w grid
1752                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
1753                                    get_topography_top_index_ji( j, i, 'w' ) )
1754             ENDDO
1755          ENDDO
1756          nzt_topo_nestbc_s = nzt_topo_nestbc_s + 1
1757       ENDIF
1758
1759       nzt_topo_nestbc_n = nzb
1760       IF ( nest_bound_n )  THEN
1761          j = nyn + 1
1762          DO  i = nxl, nxr
1763!
1764!--             Concept need to be reconsidered for 3D-topography
1765!--             Determine largest topography index on scalar grid
1766                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
1767                                    get_topography_top_index_ji( j, i, 's' ) )
1768!
1769!--             Determine largest topography index on u grid
1770                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
1771                                    get_topography_top_index_ji( j, i, 'u' ) )
1772!
1773!--             Determine largest topography index on v grid
1774                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
1775                                    get_topography_top_index_ji( j, i, 'v' ) )
1776!
1777!--             Determine largest topography index on w grid
1778                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
1779                                    get_topography_top_index_ji( j, i, 'w' ) )
1780          ENDDO
1781          nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
1782       ENDIF
1783
1784#if defined( __parallel )
1785!
1786!--       Determine global topography-top index along child boundary.
1787          dum = nzb
1788          CALL MPI_ALLREDUCE( nzt_topo_nestbc_l, dum, 1, MPI_INTEGER,          &
1789                              MPI_MAX, comm1dy, ierr )
1790          nzt_topo_nestbc_l = dum
1791
1792          dum = nzb
1793          CALL MPI_ALLREDUCE( nzt_topo_nestbc_r, dum, 1, MPI_INTEGER,          &
1794                              MPI_MAX, comm1dy, ierr )
1795          nzt_topo_nestbc_r = dum
1796
1797          dum = nzb
1798          CALL MPI_ALLREDUCE( nzt_topo_nestbc_n, dum, 1, MPI_INTEGER,          &
1799                              MPI_MAX, comm1dx, ierr )
1800          nzt_topo_nestbc_n = dum
1801
1802          dum = nzb
1803          CALL MPI_ALLREDUCE( nzt_topo_nestbc_s, dum, 1, MPI_INTEGER,          &
1804                              MPI_MAX, comm1dx, ierr )
1805          nzt_topo_nestbc_s = dum
1806#endif
1807!
1808!--    Then determine the maximum number of near-wall nodes per wall point based
1809!--    on the grid-spacing ratios.
1810       nzt_topo_max = MAX( nzt_topo_nestbc_l, nzt_topo_nestbc_r,               &
1811                           nzt_topo_nestbc_s, nzt_topo_nestbc_n )
1812!
1813!--    Note that the outer division must be integer division.
1814       ni = CEILING( cg%dx / dx ) / 2
1815       nj = CEILING( cg%dy / dy ) / 2
1816       nk = 1
1817       DO  k = 1, nzt_topo_max
1818          nk = MAX( nk, CEILING( cg%dzu(kco(k)+1) / dzu(k) ) )
1819       ENDDO
1820       nk = nk / 2   !  Note that this must be integer division.
1821       ncorr =  MAX( ni, nj, nk )
1822
1823       ALLOCATE( lcr(0:ncorr-1) )
1824       lcr = 1.0_wp
1825
1826       z0_topo = roughness_length
1827!
1828!--    First horizontal walls. Note that also logc_w_? and logc_ratio_w_? and
1829!--    logc_kbounds_* need to be allocated and initialized here.
1830!--    Left boundary
1831       IF ( nest_bound_l )  THEN
1832
1833          ALLOCATE( logc_u_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
1834          ALLOCATE( logc_v_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
1835          ALLOCATE( logc_w_l(1:2,nzb:nzt_topo_nestbc_l,nys:nyn) )
1836          ALLOCATE( logc_kbounds_u_l(1:2,nys:nyn) )
1837          ALLOCATE( logc_kbounds_v_l(1:2,nys:nyn) )
1838          ALLOCATE( logc_kbounds_w_l(1:2,nys:nyn) )
1839          ALLOCATE( logc_ratio_u_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
1840          ALLOCATE( logc_ratio_v_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
1841          ALLOCATE( logc_ratio_w_l(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_l,nys:nyn) )
1842          logc_u_l       = 0
1843          logc_v_l       = 0
1844          logc_w_l       = 0
1845          logc_ratio_u_l = 1.0_wp
1846          logc_ratio_v_l = 1.0_wp
1847          logc_ratio_w_l = 1.0_wp
1848          direction      = 1
1849          inc            = 1
1850
1851          DO  j = nys, nyn
1852!
1853!--          Left boundary for u
1854             i   = 0
1855!
1856!--          For loglaw correction the roughness z0 is required. z0, however,
1857!--          is part of the surfacetypes now. Set default roughness instead.
1858!--          Determine topography top index on u-grid
1859             kb  = get_topography_top_index_ji( j, i, 'u' )
1860             k   = kb + 1
1861             wall_index = kb
1862
1863             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
1864                            j, inc, wall_index, z0_topo,                       &
1865                            kb, direction, ncorr )
1866
1867             logc_u_l(1,k,j) = lc
1868             logc_ratio_u_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
1869             lcr(0:ncorr-1) = 1.0_wp
1870!
1871!--          Left boundary for v
1872             i   = -1
1873!
1874!--          Determine topography top index on v-grid
1875             kb  = get_topography_top_index_ji( j, i, 'v' )
1876             k   = kb + 1
1877             wall_index = kb
1878
1879             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
1880                            j, inc, wall_index, z0_topo,                       &
1881                            kb, direction, ncorr )
1882
1883             logc_v_l(1,k,j) = lc
1884             logc_ratio_v_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
1885             lcr(0:ncorr-1) = 1.0_wp
1886
1887          ENDDO
1888
1889       ENDIF
1890!
1891!--    Right boundary
1892       IF ( nest_bound_r )  THEN
1893           
1894          ALLOCATE( logc_u_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
1895          ALLOCATE( logc_v_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )
1896          ALLOCATE( logc_w_r(1:2,nzb:nzt_topo_nestbc_r,nys:nyn) )         
1897          ALLOCATE( logc_kbounds_u_r(1:2,nys:nyn) )
1898          ALLOCATE( logc_kbounds_v_r(1:2,nys:nyn) )
1899          ALLOCATE( logc_kbounds_w_r(1:2,nys:nyn) )
1900          ALLOCATE( logc_ratio_u_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
1901          ALLOCATE( logc_ratio_v_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
1902          ALLOCATE( logc_ratio_w_r(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_r,nys:nyn) )
1903          logc_u_r       = 0
1904          logc_v_r       = 0
1905          logc_w_r       = 0
1906          logc_ratio_u_r = 1.0_wp
1907          logc_ratio_v_r = 1.0_wp
1908          logc_ratio_w_r = 1.0_wp
1909          direction      = 1
1910          inc            = 1
1911
1912          DO  j = nys, nyn
1913!
1914!--          Right boundary for u
1915             i   = nxr + 1
1916!
1917!--          For loglaw correction the roughness z0 is required. z0, however,
1918!--          is part of the surfacetypes now, so call subroutine according
1919!--          to the present surface tpye.
1920!--          Determine topography top index on u-grid
1921             kb  = get_topography_top_index_ji( j, i, 'u' )
1922             k   = kb + 1
1923             wall_index = kb
1924
1925             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
1926                                                 j, inc, wall_index, z0_topo,  &
1927                                                 kb, direction, ncorr )
1928
1929             logc_u_r(1,k,j) = lc
1930             logc_ratio_u_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
1931             lcr(0:ncorr-1) = 1.0_wp
1932!
1933!--          Right boundary for v
1934             i   = nxr + 1
1935!
1936!--          Determine topography top index on v-grid
1937             kb  = get_topography_top_index_ji( j, i, 'v' )
1938             k   = kb + 1
1939             wall_index = kb
1940
1941             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
1942                                                 j, inc, wall_index, z0_topo,  &
1943                                                 kb, direction, ncorr )
1944
1945             logc_v_r(1,k,j) = lc
1946             logc_ratio_v_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
1947             lcr(0:ncorr-1) = 1.0_wp
1948
1949          ENDDO
1950
1951       ENDIF
1952!
1953!--    South boundary
1954       IF ( nest_bound_s )  THEN
1955
1956          ALLOCATE( logc_u_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1957          ALLOCATE( logc_v_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1958          ALLOCATE( logc_w_s(1:2,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1959          ALLOCATE( logc_kbounds_u_s(1:2,nxl:nxr) )
1960          ALLOCATE( logc_kbounds_v_s(1:2,nxl:nxr) )
1961          ALLOCATE( logc_kbounds_w_s(1:2,nxl:nxr) )
1962          ALLOCATE( logc_ratio_u_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1963          ALLOCATE( logc_ratio_v_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1964          ALLOCATE( logc_ratio_w_s(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_s,nxl:nxr) )
1965          logc_u_s       = 0
1966          logc_v_s       = 0
1967          logc_w_s       = 0
1968          logc_ratio_u_s = 1.0_wp
1969          logc_ratio_v_s = 1.0_wp
1970          logc_ratio_w_s = 1.0_wp
1971          direction      = 1
1972          inc            = 1
1973
1974          DO  i = nxl, nxr
1975!
1976!--          South boundary for u
1977             j   = -1
1978!
1979!--          Determine topography top index on u-grid
1980             kb  = get_topography_top_index_ji( j, i, 'u' )
1981             k   = kb + 1
1982             wall_index = kb
1983
1984             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
1985                                                 j, inc, wall_index, z0_topo,  &
1986                                                 kb, direction, ncorr )
1987
1988             logc_u_s(1,k,i) = lc
1989             logc_ratio_u_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
1990             lcr(0:ncorr-1) = 1.0_wp
1991!
1992!--          South boundary for v
1993             j   = 0
1994!
1995!--          Determine topography top index on v-grid
1996             kb  = get_topography_top_index_ji( j, i, 'v' )
1997             k   = kb + 1
1998             wall_index = kb
1999
2000             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
2001                                                 j, inc, wall_index, z0_topo,  &
2002                                                 kb, direction, ncorr )
2003
2004             logc_v_s(1,k,i) = lc
2005             logc_ratio_v_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2006             lcr(0:ncorr-1) = 1.0_wp
2007
2008          ENDDO
2009
2010       ENDIF
2011!
2012!--    North boundary
2013       IF ( nest_bound_n )  THEN
2014
2015          ALLOCATE( logc_u_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2016          ALLOCATE( logc_v_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2017          ALLOCATE( logc_w_n(1:2,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2018          ALLOCATE( logc_kbounds_u_n(1:2,nxl:nxr) )
2019          ALLOCATE( logc_kbounds_v_n(1:2,nxl:nxr) )
2020          ALLOCATE( logc_kbounds_w_n(1:2,nxl:nxr) )
2021          ALLOCATE( logc_ratio_u_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2022          ALLOCATE( logc_ratio_v_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2023          ALLOCATE( logc_ratio_w_n(1:2,0:ncorr-1,nzb:nzt_topo_nestbc_n,nxl:nxr) )
2024          logc_u_n       = 0
2025          logc_v_n       = 0
2026          logc_w_n       = 0
2027          logc_ratio_u_n = 1.0_wp
2028          logc_ratio_v_n = 1.0_wp
2029          logc_ratio_w_n = 1.0_wp
2030          direction      = 1
2031          inc            = 1
2032
2033          DO  i = nxl, nxr
2034!
2035!--          North boundary for u
2036             j   = nyn + 1
2037!
2038!--          Determine topography top index on u-grid
2039             kb  = get_topography_top_index_ji( j, i, 'u' )
2040             k   = kb + 1
2041             wall_index = kb
2042
2043             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
2044                                                 j, inc, wall_index, z0_topo,  &
2045                                                 kb, direction, ncorr )
2046
2047             logc_u_n(1,k,i) = lc
2048             logc_ratio_u_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2049             lcr(0:ncorr-1) = 1.0_wp
2050!
2051!--          North boundary for v
2052             j   = nyn + 1
2053!
2054!--          Determine topography top index on v-grid
2055             kb  = get_topography_top_index_ji( j, i, 'v' )
2056             k   = kb + 1
2057             wall_index = kb
2058
2059             CALL pmci_define_loglaw_correction_parameters( lc, lcr, k,        &
2060                                                 j, inc, wall_index, z0_topo,  &
2061                                                 kb, direction, ncorr )
2062             logc_v_n(1,k,i) = lc
2063             logc_ratio_v_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2064             lcr(0:ncorr-1) = 1.0_wp
2065
2066          ENDDO
2067
2068       ENDIF
2069!       
2070!--    Then vertical walls and corners if necessary
2071       IF ( topography /= 'flat' )  THEN
2072!
2073!--       Workaround, set z0 at vertical surfaces simply to the given roughness
2074!--       lenth, which is required to determine the logarithmic correction
2075!--       factors at the child boundaries, which are at the ghost-points.
2076!--       The surface data type for vertical surfaces, however, is not defined
2077!--       at ghost-points, so that no z0 can be retrieved at this point.
2078!--       Maybe, revise this later and define vertical surface datattype also
2079!--       at ghost-points.
2080          z0_topo = roughness_length
2081
2082          kb = 0       ! kb is not used when direction > 1       
2083!       
2084!--       Left boundary
2085          IF ( nest_bound_l )  THEN
2086             logc_kbounds_u_l(1:2,nys:nyn) = 0
2087             logc_kbounds_v_l(1:2,nys:nyn) = 0             
2088             logc_kbounds_w_l(1:2,nys:nyn) = 0
2089             
2090             direction  = 2
2091
2092             DO  j = nys, nyn
2093!
2094!--             Determine the lowest k-indices for u at j,i, j+1,i and j-1,i.
2095                i             = 0
2096                k_wall_u_ji   = get_topography_top_index_ji( j,   i, 'u' )
2097                k_wall_u_ji_p = get_topography_top_index_ji( j+1, i, 'u' )
2098                k_wall_u_ji_m = get_topography_top_index_ji( j-1, i, 'u' )
2099!
2100!--             Wall for u on the south side.
2101                IF ( ( k_wall_u_ji <  k_wall_u_ji_m ) .AND.                    &
2102                     ( k_wall_u_ji >= k_wall_u_ji_p ) )  THEN
2103                   inc        =  1
2104                   wall_index =  j
2105!
2106!--                Store the kbounds for use in pmci_interp_tril_lr.
2107                   logc_kbounds_u_l(1,j) = k_wall_u_ji + 1
2108                   logc_kbounds_u_l(2,j) = k_wall_u_ji_m
2109                   DO  k = logc_kbounds_u_l(1,j), logc_kbounds_u_l(2,j)
2110                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2111                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2112                           ncorr )
2113!
2114!--                   The direction of the wall-normal index is stored as the
2115!--                   sign of the logc-element.
2116                      logc_u_l(2,k,j) = inc * lc
2117                      logc_ratio_u_l(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2118                      lcr(0:ncorr-1) = 1.0_wp
2119                   ENDDO
2120                ENDIF
2121!
2122!--             Wall for u on the north side.
2123                IF ( ( k_wall_u_ji <  k_wall_u_ji_p ) .AND.                    &
2124                     ( k_wall_u_ji >= k_wall_u_ji_m ) )  THEN
2125                   inc        = -1
2126                   wall_index =  j + 1
2127!
2128!--                Store the kbounds for use in pmci_interp_tril_lr.                   
2129                   logc_kbounds_u_l(1,j) = k_wall_u_ji + 1
2130                   logc_kbounds_u_l(2,j) = k_wall_u_ji_p
2131                   DO  k = logc_kbounds_u_l(1,j), logc_kbounds_u_l(2,j)
2132                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2133                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2134                           ncorr )
2135!
2136!--                   The direction of the wall-normal index is stored as the
2137!--                   sign of the logc-element.
2138                      logc_u_l(2,k,j) = inc * lc
2139                      logc_ratio_u_l(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2140                      lcr(0:ncorr-1) = 1.0_wp
2141                   ENDDO
2142                ENDIF
2143!
2144!--             Determine the lowest k-indices for w at j,i, j+1,i and j-1,i.
2145                i             = -1
2146                k_wall_w_ji   = get_topography_top_index_ji( j,   i, 'w' )
2147                k_wall_w_ji_p = get_topography_top_index_ji( j+1, i, 'w' )
2148                k_wall_w_ji_m = get_topography_top_index_ji( j-1, i, 'w' )
2149!
2150!--             Wall for w on the south side.               
2151                IF ( ( k_wall_w_ji <  k_wall_w_ji_m ) .AND.                    &
2152                     ( k_wall_w_ji >= k_wall_w_ji_p ) )  THEN
2153                   inc        =  1
2154                   wall_index =  j
2155!
2156!--                Store the kbounds for use in pmci_interp_tril_lr.
2157                   logc_kbounds_w_l(1,j) = k_wall_w_ji + 1
2158                   logc_kbounds_w_l(2,j) = k_wall_w_ji_m
2159                   DO  k = logc_kbounds_w_l(1,j), logc_kbounds_w_l(2,j)
2160                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2161                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2162                           ncorr )
2163!
2164!--                   The direction of the wall-normal index is stored as the
2165!--                   sign of the logc-element.
2166                      logc_w_l(2,k,j) = inc * lc
2167                      logc_ratio_w_l(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2168                      lcr(0:ncorr-1) = 1.0_wp
2169                   ENDDO
2170                ENDIF
2171!
2172!--             Wall for w on the north side.
2173                IF ( ( k_wall_w_ji <  k_wall_w_ji_p ) .AND.                    &
2174                     ( k_wall_w_ji >= k_wall_w_ji_m ) )  THEN
2175                   inc        = -1
2176                   wall_index =  j+1
2177!
2178!--                Store the kbounds for use in pmci_interp_tril_lr.
2179                   logc_kbounds_w_l(1,j) = k_wall_w_ji + 1
2180                   logc_kbounds_w_l(2,j) = k_wall_w_ji_p
2181                   DO  k = logc_kbounds_w_l(1,j), logc_kbounds_w_l(2,j)
2182                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2183                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2184                           ncorr )
2185!
2186!--                   The direction of the wall-normal index is stored as the
2187!--                   sign of the logc-element.
2188                      logc_w_l(2,k,j) = inc * lc
2189                      logc_ratio_w_l(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2190                      lcr(0:ncorr-1) = 1.0_wp
2191                   ENDDO
2192                ENDIF
2193                   
2194             ENDDO
2195
2196          ENDIF   !  IF ( nest_bound_l )
2197!       
2198!--       Right boundary
2199          IF ( nest_bound_r )  THEN
2200             logc_kbounds_u_r(1:2,nys:nyn) = 0
2201             logc_kbounds_v_r(1:2,nys:nyn) = 0             
2202             logc_kbounds_w_r(1:2,nys:nyn) = 0
2203
2204             direction  = 2
2205             i  = nx + 1
2206
2207             DO  j = nys, nyn
2208!
2209!--             Determine the lowest k-indices for u at j,i, j+1,i and j-1,i.
2210                k_wall_u_ji   = get_topography_top_index_ji( j,   i, 'u' )
2211                k_wall_u_ji_p = get_topography_top_index_ji( j+1, i, 'u' )
2212                k_wall_u_ji_m = get_topography_top_index_ji( j-1, i, 'u' )
2213!
2214!--             Wall for u on the south side.
2215                IF ( ( k_wall_u_ji <  k_wall_u_ji_m ) .AND.                    &
2216                     ( k_wall_u_ji >= k_wall_u_ji_p ) )  THEN
2217                   inc        =  1
2218                   wall_index =  j
2219!
2220!--                Store the kbounds for use in pmci_interp_tril_lr.                 
2221                   logc_kbounds_u_r(1,j) = k_wall_u_ji + 1
2222                   logc_kbounds_u_r(2,j) = k_wall_u_ji_m
2223                   DO  k = logc_kbounds_u_r(1,j), logc_kbounds_u_r(2,j)
2224                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2225                           k, j, inc, wall_index, z0_topo, kb, direction, ncorr )
2226!
2227!--                   The direction of the wall-normal index is stored as the
2228!--                   sign of the logc-element.
2229                      logc_u_r(2,k,j) = inc * lc
2230                      logc_ratio_u_r(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2231                      lcr(0:ncorr-1) = 1.0_wp
2232                   ENDDO
2233                ENDIF
2234!
2235!--             Wall for u on the south side.
2236                IF ( ( k_wall_u_ji <  k_wall_u_ji_p ) .AND.                    &
2237                     ( k_wall_u_ji >= k_wall_u_ji_m ) )  THEN
2238                   inc        = -1
2239                   wall_index =  j + 1                 
2240!
2241!--                Store the kbounds for use in pmci_interp_tril_lr.                   
2242                   logc_kbounds_u_r(1,j) = k_wall_u_ji + 1
2243                   logc_kbounds_u_r(2,j) = k_wall_u_ji_p
2244                   DO  k = logc_kbounds_u_r(1,j), logc_kbounds_u_r(2,j)
2245                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2246                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2247                           ncorr )
2248!
2249!--                   The direction of the wall-normal index is stored as the
2250!--                   sign of the logc-element.
2251                      logc_u_r(2,k,j) = inc * lc
2252                      logc_ratio_u_r(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2253                      lcr(0:ncorr-1) = 1.0_wp
2254                   ENDDO
2255                ENDIF
2256!
2257!--             Determine the lowest k-indices for w at j,i, j+1,i and j-1,i.
2258                k_wall_w_ji   = get_topography_top_index_ji( j,   i, 'w' )
2259                k_wall_w_ji_p = get_topography_top_index_ji( j+1, i, 'w' )
2260                k_wall_w_ji_m = get_topography_top_index_ji( j-1, i, 'w' )
2261!
2262!--             Wall for w on the south side.               
2263                IF ( ( k_wall_w_ji <  k_wall_w_ji_m ) .AND.                    &
2264                     ( k_wall_w_ji >= k_wall_w_ji_p ) )  THEN
2265                   inc        =  1
2266                   wall_index =  j
2267!
2268!--                Store the kbounds for use in pmci_interp_tril_lr.                   
2269                   logc_kbounds_w_r(1,j) = k_wall_w_ji + 1
2270                   logc_kbounds_w_r(2,j) = k_wall_w_ji_m
2271                   DO  k = logc_kbounds_w_r(1,j), logc_kbounds_w_r(2,j)
2272                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2273                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2274                           ncorr )
2275!
2276!--                   The direction of the wall-normal index is stored as the
2277!--                   sign of the logc-element.
2278                      logc_w_r(2,k,j) = inc * lc
2279                      logc_ratio_w_r(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2280                      lcr(0:ncorr-1) = 1.0_wp
2281                   ENDDO
2282                ENDIF
2283!
2284!--             Wall for w on the north side.
2285                IF ( ( k_wall_w_ji <  k_wall_w_ji_p ) .AND.                    &
2286                     ( k_wall_w_ji >= k_wall_w_ji_m ) )  THEN
2287                   inc        = -1
2288                   wall_index =  j+1
2289!
2290!--                Store the kbounds for use in pmci_interp_tril_lr.                   
2291                   logc_kbounds_w_r(1,j) = k_wall_w_ji + 1
2292                   logc_kbounds_w_r(2,j) = k_wall_w_ji_p
2293                   DO  k = logc_kbounds_w_r(1,j), logc_kbounds_w_r(2,j)
2294                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2295                           k, j, inc, wall_index, z0_topo, kb, direction,      &
2296                           ncorr )
2297!
2298!--                   The direction of the wall-normal index is stored as the
2299!--                   sign of the logc-element.
2300                      logc_w_r(2,k,j) = inc * lc
2301                      logc_ratio_w_r(2,0:ncorr-1,k,j) = lcr(0:ncorr-1)
2302                      lcr(0:ncorr-1) = 1.0_wp
2303                   ENDDO
2304                ENDIF
2305                   
2306             ENDDO
2307             
2308          ENDIF   !  IF ( nest_bound_r )
2309!       
2310!--       South boundary
2311          IF ( nest_bound_s )  THEN
2312             logc_kbounds_u_s(1:2,nxl:nxr) = 0
2313             logc_kbounds_v_s(1:2,nxl:nxr) = 0
2314             logc_kbounds_w_s(1:2,nxl:nxr) = 0
2315
2316             direction  = 3
2317
2318             DO  i = nxl, nxr
2319!
2320!--             Determine the lowest k-indices for v at j,i, j,i+1 and j,i-1.
2321                j             = 0               
2322                k_wall_v_ji   = get_topography_top_index_ji( j, i,   'v' )
2323                k_wall_v_ji_p = get_topography_top_index_ji( j, i+1, 'v' )
2324                k_wall_v_ji_m = get_topography_top_index_ji( j, i-1, 'v' )
2325!
2326!--             Wall for v on the left side.
2327                IF ( ( k_wall_v_ji <  k_wall_v_ji_m ) .AND.                    &
2328                     ( k_wall_v_ji >= k_wall_v_ji_p ) )  THEN
2329                   inc        =  1
2330                   wall_index =  i
2331!
2332!--                Store the kbounds for use in pmci_interp_tril_sn.                   
2333                   logc_kbounds_v_s(1,i) = k_wall_v_ji + 1
2334                   logc_kbounds_v_s(2,i) = k_wall_v_ji_m
2335                   DO  k = logc_kbounds_v_s(1,i), logc_kbounds_v_s(2,i)
2336                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2337                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2338                           ncorr )
2339!
2340!--                   The direction of the wall-normal index is stored as the
2341!--                   sign of the logc-element.
2342                      logc_v_s(2,k,i) = inc * lc
2343                      logc_ratio_v_s(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2344                      lcr(0:ncorr-1) = 1.0_wp
2345                   ENDDO
2346                ENDIF
2347!
2348!--             Wall for v on the right side.
2349                IF ( ( k_wall_v_ji <  k_wall_v_ji_p ) .AND.                    &
2350                     ( k_wall_v_ji >= k_wall_v_ji_m ) )  THEN
2351                   inc        = -1
2352                   wall_index =  i+1
2353!
2354!--                Store the kbounds for use in pmci_interp_tril_sn.                   
2355                   logc_kbounds_v_s(1,i) = k_wall_v_ji + 1
2356                   logc_kbounds_v_s(2,i) = k_wall_v_ji_p
2357                   DO  k = logc_kbounds_v_s(1,i), logc_kbounds_v_s(2,i)
2358                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2359                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2360                           ncorr )
2361!
2362!--                   The direction of the wall-normal index is stored as the
2363!--                   sign of the logc-element.
2364                      logc_v_s(2,k,i) = inc * lc
2365                      logc_ratio_v_s(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2366                      lcr(0:ncorr-1) = 1.0_wp
2367                   ENDDO
2368                ENDIF
2369!
2370!--             Determine the lowest k-indices for w at j,i, j,i+1 and j,i-1.
2371                j             = -1
2372                k_wall_w_ji   = get_topography_top_index_ji( j, i,   'w' )
2373                k_wall_w_ji_p = get_topography_top_index_ji( j, i+1, 'w' )
2374                k_wall_w_ji_m = get_topography_top_index_ji( j, i-1, 'w' )
2375!
2376!--             Wall for w on the left side.
2377                IF ( ( k_wall_w_ji <  k_wall_w_ji_m ) .AND.                    &
2378                     ( k_wall_w_ji >= k_wall_w_ji_p ) )  THEN
2379                   inc        =  1
2380                   wall_index =  i
2381!
2382!--                Store the kbounds for use in pmci_interp_tril_sn.
2383                   logc_kbounds_w_s(1,i) = k_wall_w_ji + 1
2384                   logc_kbounds_w_s(2,i) = k_wall_w_ji_m
2385                   DO  k = logc_kbounds_w_s(1,i), logc_kbounds_w_s(2,i)
2386                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2387                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2388                           ncorr )
2389!
2390!--                   The direction of the wall-normal index is stored as the
2391!--                   sign of the logc-element.
2392                      logc_w_s(2,k,i) = inc * lc
2393                      logc_ratio_w_s(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2394                      lcr(0:ncorr-1) = 1.0_wp
2395                   ENDDO
2396                ENDIF
2397!
2398!--             Wall for w on the right side.
2399                IF ( ( k_wall_w_ji <  k_wall_w_ji_p ) .AND.                    &
2400                     ( k_wall_w_ji >= k_wall_w_ji_m ) )  THEN
2401                   inc        = -1
2402                   wall_index =  i+1
2403!
2404!--                Store the kbounds for use in pmci_interp_tril_sn.
2405                   logc_kbounds_w_s(1,i) = k_wall_w_ji + 1
2406                   logc_kbounds_w_s(2,i) = k_wall_w_ji_p
2407                   DO  k = logc_kbounds_w_s(1,i), logc_kbounds_w_s(2,i)
2408                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2409                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2410                           ncorr )
2411!
2412!--                   The direction of the wall-normal index is stored as the
2413!--                   sign of the logc-element.
2414                      logc_w_s(2,k,i) = inc * lc
2415                      logc_ratio_w_s(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2416                      lcr(0:ncorr-1) = 1.0_wp
2417                   ENDDO
2418                ENDIF
2419
2420             ENDDO
2421
2422          ENDIF   !  IF (nest_bound_s )
2423!       
2424!--       North boundary
2425          IF ( nest_bound_n )  THEN
2426             logc_kbounds_u_n(1:2,nxl:nxr) = 0             
2427             logc_kbounds_v_n(1:2,nxl:nxr) = 0
2428             logc_kbounds_w_n(1:2,nxl:nxr) = 0
2429
2430             direction  = 3
2431             j  = ny + 1
2432
2433             DO  i = nxl, nxr
2434!
2435!--             Determine the lowest k-indices for v at j,i, j,i+1 and j,i-1
2436                k_wall_v_ji   = get_topography_top_index_ji( j, i,   'v' )
2437                k_wall_v_ji_p = get_topography_top_index_ji( j, i+1, 'v' )
2438                k_wall_v_ji_m = get_topography_top_index_ji( j, i-1, 'v' )
2439!
2440!--             Wall for v on the left side.
2441                IF ( ( k_wall_v_ji <  k_wall_v_ji_m ) .AND.                    &
2442                     ( k_wall_v_ji >= k_wall_v_ji_p ) )  THEN
2443                   inc        = 1
2444                   wall_index = i                   
2445!
2446!--                Store the kbounds for use in pmci_interp_tril_sn.
2447                   logc_kbounds_v_n(1,i) = k_wall_v_ji + 1
2448                   logc_kbounds_v_n(2,i) = k_wall_v_ji_m
2449                   DO  k = logc_kbounds_v_n(1,i), logc_kbounds_v_n(2,i)
2450                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2451                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2452                           ncorr )
2453!
2454!--                   The direction of the wall-normal index is stored as the
2455!--                   sign of the logc-element.
2456                      logc_v_n(2,k,i) = inc * lc
2457                      logc_ratio_v_n(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2458                      lcr(0:ncorr-1) = 1.0_wp
2459                   ENDDO
2460                ENDIF
2461!
2462!--             Wall for v on the right side.
2463                IF ( ( k_wall_v_ji <  k_wall_v_ji_p ) .AND.                    &
2464                     ( k_wall_v_ji >= k_wall_v_ji_m ) )  THEN
2465                   inc        = -1
2466                   wall_index =  i + 1
2467!
2468!--                Store the kbounds for use in pmci_interp_tril_sn.
2469                   logc_kbounds_v_n(1,i) = k_wall_v_ji + 1
2470                   logc_kbounds_v_n(2,i) = k_wall_v_ji_p
2471                   DO  k = logc_kbounds_v_n(1,i), logc_kbounds_v_n(2,i)
2472                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2473                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2474                           ncorr )
2475!
2476!--                   The direction of the wall-normal index is stored as the
2477!--                   sign of the logc-element.
2478                      logc_v_n(2,k,i) = inc * lc
2479                      logc_ratio_v_n(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2480                      lcr(0:ncorr-1) = 1.0_wp
2481                   ENDDO
2482                ENDIF
2483!
2484!--             Determine the lowest k-indices for w at j,i, j,i+1 and j,i-1.
2485                k_wall_w_ji   = get_topography_top_index_ji( j, i,   'w' )
2486                k_wall_w_ji_p = get_topography_top_index_ji( j, i+1, 'w' )
2487                k_wall_w_ji_m = get_topography_top_index_ji( j, i-1, 'w' )                   
2488!
2489!--             Wall for w on the left side.
2490                IF ( ( k_wall_w_ji <  k_wall_w_ji_m ) .AND.                    &
2491                     ( k_wall_w_ji >= k_wall_w_ji_p ) )  THEN
2492                   inc        = 1
2493                   wall_index = i
2494!
2495!--                Store the kbounds for use in pmci_interp_tril_sn.
2496                   logc_kbounds_w_n(1,i) = k_wall_w_ji + 1
2497                   logc_kbounds_w_n(2,i) = k_wall_w_ji_m
2498                   DO  k = logc_kbounds_w_n(1,i), logc_kbounds_w_n(2,i)
2499                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2500                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2501                           ncorr )
2502!
2503!--                   The direction of the wall-normal index is stored as the
2504!--                   sign of the logc-element.
2505                      logc_w_n(2,k,i) = inc * lc
2506                      logc_ratio_w_n(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2507                      lcr(0:ncorr-1) = 1.0_wp
2508                   ENDDO
2509                ENDIF
2510!
2511!--             Wall for w on the right side, but not on the left side
2512                IF ( ( k_wall_w_ji <  k_wall_w_ji_p ) .AND.                    &
2513                     ( k_wall_w_ji >= k_wall_w_ji_m ) )  THEN
2514                   inc        = -1
2515                   wall_index =  i+1
2516!
2517!--                Store the kbounds for use in pmci_interp_tril_sn.
2518                   logc_kbounds_w_n(1,i) = k_wall_w_ji + 1
2519                   logc_kbounds_w_n(2,i) = k_wall_w_ji_p
2520                   DO  k = logc_kbounds_w_n(1,i), logc_kbounds_w_n(2,i)
2521                      CALL pmci_define_loglaw_correction_parameters( lc, lcr,  &
2522                           k, i, inc, wall_index, z0_topo, kb, direction,      &
2523                           ncorr )
2524!
2525!--                   The direction of the wall-normal index is stored as the
2526!--                   sign of the logc-element.
2527                      logc_w_n(2,k,i) = inc * lc
2528                      logc_ratio_w_n(2,0:ncorr-1,k,i) = lcr(0:ncorr-1)
2529                      lcr(0:ncorr-1) = 1.0_wp
2530                   ENDDO
2531                ENDIF
2532
2533             ENDDO
2534
2535          ENDIF   !  IF ( nest_bound_n )
2536
2537       ENDIF   !  IF ( topography /= 'flat' )
2538
2539    END SUBROUTINE pmci_init_loglaw_correction
2540
2541
2542
2543    SUBROUTINE pmci_define_loglaw_correction_parameters( lc, lcr, k, ij, inc,  &
2544         wall_index, z0_l, kb, direction, ncorr )
2545
2546       IMPLICIT NONE
2547
2548       INTEGER(iwp), INTENT(IN)  ::  direction  !<
2549       INTEGER(iwp), INTENT(IN)  ::  ij         !<
2550       INTEGER(iwp), INTENT(IN)  ::  inc        !<
2551       INTEGER(iwp), INTENT(IN)  ::  k          !<
2552       INTEGER(iwp), INTENT(IN)  ::  kb         !<
2553       INTEGER(iwp), INTENT(OUT) ::  lc         !<
2554       INTEGER(iwp), INTENT(IN)  ::  ncorr      !<
2555       INTEGER(iwp), INTENT(IN)  ::  wall_index !<
2556
2557       INTEGER(iwp) ::  alcorr       !<
2558       INTEGER(iwp) ::  corr_index   !<
2559       INTEGER(iwp) ::  lcorr        !<
2560
2561       LOGICAL      ::  more         !<
2562
2563       REAL(wp), DIMENSION(0:ncorr-1), INTENT(OUT) ::  lcr     !<
2564       REAL(wp), INTENT(IN)      ::  z0_l                      !<
2565     
2566       REAL(wp)     ::  logvelc1     !<
2567     
2568
2569       SELECT CASE ( direction )
2570
2571          CASE (1)   !  k
2572             more = .TRUE.
2573             lcorr = 0
2574             DO  WHILE ( more .AND. lcorr <= ncorr-1 )
2575                corr_index = k + lcorr
2576                IF ( lcorr == 0 )  THEN
2577                   CALL pmci_find_logc_pivot_k( lc, logvelc1, z0_l, kb )
2578                ENDIF
2579               
2580                IF ( corr_index < lc )  THEN
2581                   lcr(lcorr) = LOG( ( zu(k) - zw(kb) ) / z0_l ) / logvelc1
2582                   more = .TRUE.
2583                ELSE
2584                   lcr(lcorr) = 1.0
2585                   more = .FALSE.
2586                ENDIF
2587                lcorr = lcorr + 1
2588             ENDDO
2589
2590          CASE (2)   !  j
2591             more = .TRUE.
2592             lcorr  = 0
2593             alcorr = 0
2594             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
2595                corr_index = ij + lcorr   ! In this case (direction = 2) ij is j
2596                IF ( lcorr == 0 )  THEN
2597                   CALL pmci_find_logc_pivot_j( lc, logvelc1, ij, wall_index,  &
2598                                                z0_l, inc )
2599                ENDIF
2600!
2601!--             The role of inc here is to make the comparison operation "<"
2602!--             valid in both directions
2603                IF ( inc * corr_index < inc * lc )  THEN
2604                   lcr(alcorr) = LOG( ABS( coord_y(corr_index) + 0.5_wp * dy   &
2605                                         - coord_y(wall_index) ) / z0_l )      &
2606                                 / logvelc1
2607                   more = .TRUE.
2608                ELSE
2609                   lcr(alcorr) = 1.0_wp
2610                   more = .FALSE.
2611                ENDIF
2612                lcorr  = lcorr + inc
2613                alcorr = ABS( lcorr )
2614             ENDDO
2615
2616          CASE (3)   !  i
2617             more = .TRUE.
2618             lcorr  = 0
2619             alcorr = 0
2620             DO  WHILE ( more  .AND.  alcorr <= ncorr-1 )
2621                corr_index = ij + lcorr   ! In this case (direction = 3) ij is i
2622                IF ( lcorr == 0 )  THEN
2623                   CALL pmci_find_logc_pivot_i( lc, logvelc1, ij, wall_index,  &
2624                                                z0_l, inc )
2625                ENDIF
2626!
2627!--             The role of inc here is to make the comparison operation "<"
2628!--             valid in both directions
2629                IF ( inc * corr_index < inc * lc )  THEN
2630                   lcr(alcorr) = LOG( ABS( coord_x(corr_index) + 0.5_wp * dx   &
2631                                         - coord_x(wall_index) ) / z0_l )      &
2632                                 / logvelc1
2633                   more = .TRUE.
2634                ELSE
2635                   lcr(alcorr) = 1.0_wp
2636                   more = .FALSE.
2637                ENDIF
2638                lcorr  = lcorr + inc
2639                alcorr = ABS( lcorr )
2640             ENDDO
2641
2642       END SELECT
2643
2644    END SUBROUTINE pmci_define_loglaw_correction_parameters
2645
2646
2647
2648    SUBROUTINE pmci_find_logc_pivot_k( lc, logzc1, z0_l, kb )
2649!
2650!--    Finds the pivot node and the log-law factor for near-wall nodes for
2651!--    which the wall-parallel velocity components will be log-law corrected
2652!--    after interpolation. This subroutine is only for horizontal walls.
2653
2654       IMPLICIT NONE
2655
2656       INTEGER(iwp), INTENT(IN)  ::  kb   !<
2657       INTEGER(iwp), INTENT(OUT) ::  lc   !<
2658
2659       INTEGER(iwp) ::  kbc    !<
2660       INTEGER(iwp) ::  k1     !<
2661
2662       REAL(wp), INTENT(OUT) ::  logzc1     !<
2663       REAL(wp), INTENT(IN) ::  z0_l       !<
2664
2665       REAL(wp) ::  zuc1   !<
2666
2667
2668       kbc = nzb + 1
2669!
2670!--    kbc is the first coarse-grid point above the surface
2671       DO  WHILE ( cg%zu(kbc) < zu(kb) )
2672          kbc = kbc + 1
2673       ENDDO
2674       zuc1  = cg%zu(kbc)
2675       k1    = kb + 1
2676       DO  WHILE ( zu(k1) < zuc1 )  !  Important: must be <, not <=
2677          k1 = k1 + 1
2678       ENDDO
2679       logzc1 = LOG( (zu(k1) - zw(kb) ) / z0_l )
2680       lc = k1
2681
2682    END SUBROUTINE pmci_find_logc_pivot_k
2683
2684
2685
2686    SUBROUTINE pmci_find_logc_pivot_j( lc, logyc1, j, jw, z0_l, inc )
2687!
2688!--    Finds the pivot node and te log-law factor for near-wall nodes for
2689!--    which the wall-parallel velocity components will be log-law corrected
2690!--    after interpolation. This subroutine is only for vertical walls on
2691!--    south/north sides of the node.
2692
2693       IMPLICIT NONE
2694
2695       INTEGER(iwp), INTENT(IN)  ::  inc    !<  increment must be 1 or -1.
2696       INTEGER(iwp), INTENT(IN)  ::  j      !<
2697       INTEGER(iwp), INTENT(IN)  ::  jw     !<
2698       INTEGER(iwp), INTENT(OUT) ::  lc     !<
2699
2700       INTEGER(iwp) ::  j1       !<
2701
2702       REAL(wp), INTENT(IN) ::  z0_l   !<
2703
2704       REAL(wp) ::  logyc1   !<
2705       REAL(wp) ::  yc1      !<
2706       
2707!
2708!--    yc1 is the y-coordinate of the first coarse-grid u- and w-nodes out from
2709!--    the wall
2710       yc1  = coord_y(jw) + 0.5_wp * inc * cg%dy
2711!
2712!--    j1 is the first fine-grid index further away from the wall than yc1
2713       j1 = j
2714!
2715!--    Important: must be <, not <=
2716       DO  WHILE ( inc * ( coord_y(j1) + 0.5_wp * dy ) < inc * yc1 )
2717          j1 = j1 + inc
2718       ENDDO
2719
2720       logyc1 = LOG( ABS( coord_y(j1) + 0.5_wp * dy - coord_y(jw) ) / z0_l )
2721       lc = j1
2722
2723    END SUBROUTINE pmci_find_logc_pivot_j
2724
2725
2726
2727    SUBROUTINE pmci_find_logc_pivot_i( lc, logxc1, i, iw, z0_l, inc )
2728!
2729!--    Finds the pivot node and the log-law factor for near-wall nodes for
2730!--    which the wall-parallel velocity components will be log-law corrected
2731!--    after interpolation. This subroutine is only for vertical walls on
2732!--    south/north sides of the node.
2733
2734       IMPLICIT NONE
2735
2736       INTEGER(iwp), INTENT(IN)  ::  i      !<
2737       INTEGER(iwp), INTENT(IN)  ::  inc    !< increment must be 1 or -1.
2738       INTEGER(iwp), INTENT(IN)  ::  iw     !<
2739       INTEGER(iwp), INTENT(OUT) ::  lc     !<
2740
2741       INTEGER(iwp) ::  i1       !<
2742
2743       REAL(wp), INTENT(IN) ::  z0_l   !<
2744
2745       REAL(wp) ::  logxc1   !<
2746       REAL(wp) ::  xc1      !<
2747
2748!
2749!--    xc1 is the x-coordinate of the first coarse-grid v- and w-nodes out from
2750!--    the wall
2751       xc1  = coord_x(iw) + 0.5_wp * inc * cg%dx
2752!
2753!--    i1 is the first fine-grid index futher away from the wall than xc1.
2754       i1 = i
2755!
2756!--    Important: must be <, not <=
2757       DO  WHILE ( inc * ( coord_x(i1) + 0.5_wp *dx ) < inc * xc1 )
2758          i1 = i1 + inc
2759       ENDDO
2760     
2761       logxc1 = LOG( ABS( coord_x(i1) + 0.5_wp*dx - coord_x(iw) ) / z0_l )
2762       lc = i1
2763
2764    END SUBROUTINE pmci_find_logc_pivot_i
2765
2766
2767
2768
2769    SUBROUTINE pmci_init_anterp_tophat
2770!
2771!--    Precomputation of the child-array indices for
2772!--    corresponding coarse-grid array index and the
2773!--    Under-relaxation coefficients to be used by anterp_tophat.
2774
2775       IMPLICIT NONE
2776
2777       INTEGER(iwp) ::  i        !< Fine-grid index
2778       INTEGER(iwp) ::  ifc_o    !<
2779       INTEGER(iwp) ::  ifc_u    !<
2780       INTEGER(iwp) ::  ii       !< Coarse-grid index
2781       INTEGER(iwp) ::  istart   !<
2782       INTEGER(iwp) ::  ir       !<
2783       INTEGER(iwp) ::  j        !< Fine-grid index
2784       INTEGER(iwp) ::  jj       !< Coarse-grid index
2785       INTEGER(iwp) ::  jstart   !<
2786       INTEGER(iwp) ::  jr       !<
2787       INTEGER(iwp) ::  k        !< Fine-grid index
2788       INTEGER(iwp) ::  kk       !< Coarse-grid index
2789       INTEGER(iwp) ::  kstart   !<
2790       REAL(wp)     ::  xi       !<
2791       REAL(wp)     ::  eta      !<
2792       REAL(wp)     ::  zeta     !<
2793     
2794!
2795!--    Default values for under-relaxation lengths:
2796       IF ( anterp_relax_length_l < 0.0_wp )  THEN
2797          anterp_relax_length_l = 0.1_wp * ( nx + 1 ) * dx
2798       ENDIF
2799       IF ( anterp_relax_length_r < 0.0_wp )  THEN
2800          anterp_relax_length_r = 0.1_wp * ( nx + 1 ) * dx
2801       ENDIF
2802       IF ( anterp_relax_length_s < 0.0_wp )  THEN
2803          anterp_relax_length_s = 0.1_wp * ( ny + 1 ) * dy
2804       ENDIF
2805       IF ( anterp_relax_length_n < 0.0_wp )  THEN
2806          anterp_relax_length_n = 0.1_wp * ( ny + 1 ) * dy
2807       ENDIF
2808       IF ( anterp_relax_length_t < 0.0_wp )  THEN
2809          anterp_relax_length_t = 0.1_wp * zu(nzt)
2810       ENDIF
2811!
2812!--    First determine kctu and kctw that are the coarse-grid upper bounds for
2813!--    index k
2814       kk = 0
2815       DO  WHILE ( cg%zu(kk) <= zu(nzt) )
2816          kk = kk + 1
2817       ENDDO
2818       kctu = kk - 1
2819
2820       kk = 0
2821       DO  WHILE ( cg%zw(kk) <= zw(nzt-1) )
2822          kk = kk + 1
2823       ENDDO
2824       kctw = kk - 1
2825
2826       ALLOCATE( iflu(icl:icr) )
2827       ALLOCATE( iflo(icl:icr) )
2828       ALLOCATE( ifuu(icl:icr) )
2829       ALLOCATE( ifuo(icl:icr) )
2830       ALLOCATE( jflv(jcs:jcn) )
2831       ALLOCATE( jflo(jcs:jcn) )
2832       ALLOCATE( jfuv(jcs:jcn) )
2833       ALLOCATE( jfuo(jcs:jcn) )
2834       ALLOCATE( kflw(0:kctw) )
2835       ALLOCATE( kflo(0:kctu) )
2836       ALLOCATE( kfuw(0:kctw) )
2837       ALLOCATE( kfuo(0:kctu) )
2838
2839       ALLOCATE( ijfc_u(jcs:jcn,icl:icr) )
2840       ALLOCATE( ijfc_v(jcs:jcn,icl:icr) )
2841       ALLOCATE( ijfc_s(jcs:jcn,icl:icr) )
2842       ALLOCATE( kfc_w(0:kctw) )
2843       ALLOCATE( kfc_s(0:kctu) )
2844!
2845!--    i-indices of u for each ii-index value
2846!--    ii=icr is redundant for anterpolation
2847       istart = nxlg
2848       DO  ii = icl, icr-1
2849          i = istart
2850          DO  WHILE ( ( coord_x(i) < cg%coord_x(ii) - 0.5_wp * cg%dx )  .AND.  &
2851                      ( i < nxrg ) )
2852             i  = i + 1
2853          ENDDO
2854          iflu(ii) = MIN( MAX( i, nxlg ), nxrg )
2855          ir = i
2856          DO  WHILE ( ( coord_x(ir) <= cg%coord_x(ii) + 0.5_wp * cg%dx )  .AND.&
2857                      ( i < nxrg+1 ) )
2858             i  = i + 1
2859             ir = MIN( i, nxrg )
2860          ENDDO
2861          ifuu(ii) = MIN( MAX( i-1, iflu(ii) ), nxrg )
2862          istart = iflu(ii)
2863       ENDDO
2864       iflu(icr) = nxrg
2865       ifuu(icr) = nxrg
2866!
2867!--    i-indices of others for each ii-index value
2868!--    ii=icr is redundant for anterpolation
2869       istart = nxlg
2870       DO  ii = icl, icr-1
2871          i = istart
2872          DO  WHILE ( ( coord_x(i) + 0.5_wp * dx < cg%coord_x(ii) )  .AND.     &
2873                      ( i < nxrg ) )
2874             i  = i + 1
2875          ENDDO
2876          iflo(ii) = MIN( MAX( i, nxlg ), nxrg )
2877          ir = i
2878          DO  WHILE ( ( coord_x(ir) + 0.5_wp * dx <= cg%coord_x(ii) + cg%dx )  &
2879                      .AND.  ( i < nxrg+1 ) )
2880             i  = i + 1
2881             ir = MIN( i, nxrg )
2882          ENDDO
2883          ifuo(ii) = MIN( MAX( i-1, iflo(ii) ), nxrg )
2884          istart = iflo(ii)
2885       ENDDO
2886       iflo(icr) = nxrg
2887       ifuo(icr) = nxrg
2888!
2889!--    j-indices of v for each jj-index value
2890!--    jj=jcn is redundant for anterpolation
2891       jstart = nysg
2892       DO  jj = jcs, jcn-1
2893          j = jstart
2894          DO  WHILE ( ( coord_y(j) < cg%coord_y(jj) - 0.5_wp * cg%dy )  .AND.  &
2895                      ( j < nyng ) )
2896             j  = j + 1
2897          ENDDO
2898          jflv(jj) = MIN( MAX( j, nysg ), nyng )
2899          jr = j
2900          DO  WHILE ( ( coord_y(jr) <= cg%coord_y(jj) + 0.5_wp * cg%dy )  .AND.&
2901                      ( j < nyng+1 ) )
2902             j  = j + 1
2903             jr = MIN( j, nyng )
2904          ENDDO
2905          jfuv(jj) = MIN( MAX( j-1, jflv(jj) ), nyng )
2906          jstart = jflv(jj)
2907       ENDDO
2908       jflv(jcn) = nyng
2909       jfuv(jcn) = nyng
2910!
2911!--    j-indices of others for each jj-index value
2912!--    jj=jcn is redundant for anterpolation
2913       jstart = nysg
2914       DO  jj = jcs, jcn-1
2915          j = jstart
2916          DO  WHILE ( ( coord_y(j) + 0.5_wp * dy < cg%coord_y(jj) )  .AND.     &
2917                      ( j < nyng ) )
2918             j  = j + 1
2919          ENDDO
2920          jflo(jj) = MIN( MAX( j, nysg ), nyng )
2921          jr = j
2922          DO  WHILE ( ( coord_y(jr) + 0.5_wp * dy <= cg%coord_y(jj) + cg%dy )  &
2923                      .AND.  ( j < nyng+1 ) )
2924             j  = j + 1
2925             jr = MIN( j, nyng )
2926          ENDDO
2927          jfuo(jj) = MIN( MAX( j-1, jflo(jj) ), nyng )
2928          jstart = jflo(jj)
2929       ENDDO
2930       jflo(jcn) = nyng
2931       jfuo(jcn) = nyng
2932!
2933!--    k-indices of w for each kk-index value
2934       kstart  = 0
2935       kflw(0) = 0
2936       kfuw(0) = 0
2937       DO  kk = 1, kctw
2938          k = kstart
2939          DO  WHILE ( ( zw(k) < cg%zu(kk) )  .AND.  ( k < nzt ) )
2940             k = k + 1
2941          ENDDO
2942          kflw(kk) = MIN( MAX( k, 1 ), nzt + 1 )
2943          DO  WHILE ( ( zw(k) <= cg%zu(kk+1) )  .AND.  ( k < nzt+1 ) )
2944             k  = k + 1
2945          ENDDO
2946          kfuw(kk) = MIN( MAX( k-1, kflw(kk) ), nzt + 1 )
2947          kstart = kflw(kk)
2948       ENDDO
2949!
2950!--    k-indices of others for each kk-index value
2951       kstart  = 0
2952       kflo(0) = 0
2953       kfuo(0) = 0
2954       DO  kk = 1, kctu
2955          k = kstart
2956          DO  WHILE ( ( zu(k) < cg%zw(kk-1) )  .AND.  ( k < nzt ) )
2957             k = k + 1
2958          ENDDO
2959          kflo(kk) = MIN( MAX( k, 1 ), nzt + 1 )
2960          DO  WHILE ( ( zu(k) <= cg%zw(kk) )  .AND.  ( k < nzt+1 ) )
2961             k = k + 1
2962          ENDDO
2963          kfuo(kk) = MIN( MAX( k-1, kflo(kk) ), nzt + 1 )
2964          kstart = kflo(kk)
2965       ENDDO
2966!
2967!--    Precomputation of number of fine-grid nodes inside coarse-grid ij-faces.
2968!--    Note that ii, jj, and kk are coarse-grid indices.
2969!--    This information is needed in anterpolation.
2970       DO  ii = icl, icr
2971          ifc_u = ifuu(ii) - iflu(ii) + 1
2972          ifc_o = ifuo(ii) - iflo(ii) + 1
2973          DO  jj = jcs, jcn
2974             ijfc_u(jj,ii) = ifc_u * ( jfuo(jj) - jflo(jj) + 1 )
2975             ijfc_v(jj,ii) = ifc_o * ( jfuv(jj) - jflv(jj) + 1 )
2976             ijfc_s(jj,ii) = ifc_o * ( jfuo(jj) - jflo(jj) + 1 )             
2977          ENDDO
2978       ENDDO
2979       DO kk = 0, kctw
2980          kfc_w(kk) = kfuw(kk) - kflw(kk) + 1
2981       ENDDO
2982       DO kk = 0, kctu
2983          kfc_s(kk) = kfuo(kk) - kflo(kk) + 1
2984       ENDDO
2985!
2986!--    Spatial under-relaxation coefficients
2987       ALLOCATE( frax(icl:icr) )
2988       ALLOCATE( fray(jcs:jcn) )
2989       
2990       frax(icl:icr) = 1.0_wp
2991       fray(jcs:jcn) = 1.0_wp
2992
2993       IF ( nesting_mode /= 'vertical' )  THEN
2994          DO  ii = icl, icr
2995             IF ( ifuu(ii) < ( nx + 1 ) / 2 )  THEN   
2996                xi = ( MAX( 0.0_wp, ( cg%coord_x(ii) -                         &
2997                     lower_left_coord_x ) ) / anterp_relax_length_l )**4
2998                frax(ii) = xi / ( 1.0_wp + xi )
2999             ELSE
3000                xi = ( MAX( 0.0_wp, ( lower_left_coord_x + ( nx + 1 ) * dx -   &
3001                                      cg%coord_x(ii) ) ) /                     &
3002                       anterp_relax_length_r )**4
3003                frax(ii) = xi / ( 1.0_wp + xi )               
3004             ENDIF
3005          ENDDO
3006
3007
3008          DO  jj = jcs, jcn
3009             IF ( jfuv(jj) < ( ny + 1 ) / 2 )  THEN
3010                eta = ( MAX( 0.0_wp, ( cg%coord_y(jj) -                        &
3011                     lower_left_coord_y ) ) / anterp_relax_length_s )**4
3012                fray(jj) = eta / ( 1.0_wp + eta )
3013             ELSE
3014                eta = ( MAX( 0.0_wp, ( lower_left_coord_y + ( ny + 1 ) * dy -  &
3015                                       cg%coord_y(jj)) ) /                     &
3016                        anterp_relax_length_n )**4
3017                fray(jj) = eta / ( 1.0_wp + eta )
3018             ENDIF
3019          ENDDO
3020       ENDIF
3021     
3022       ALLOCATE( fraz(0:kctu) )
3023       DO  kk = 0, kctu
3024          zeta = ( ( zu(nzt) - cg%zu(kk) ) / anterp_relax_length_t )**4
3025          fraz(kk) = zeta / ( 1.0_wp + zeta )
3026       ENDDO
3027
3028    END SUBROUTINE pmci_init_anterp_tophat
3029
3030
3031
3032    SUBROUTINE pmci_init_tkefactor
3033
3034!
3035!--    Computes the scaling factor for the SGS TKE from coarse grid to be used
3036!--    as BC for the fine grid. Based on the Kolmogorov energy spectrum
3037!--    for the inertial subrange and assumption of sharp cut-off of the resolved
3038!--    energy spectrum. Near the surface, the reduction of TKE is made
3039!--    smaller than further away from the surface.
3040!--    Please note, in case parent and child model operate in RANS mode,
3041!--    TKE is not grid depenedent and weighting factor is one.
3042
3043       IMPLICIT NONE
3044
3045       INTEGER(iwp)        ::  k                     !< index variable along z
3046       INTEGER(iwp)        ::  k_wall                !< topography-top index along z
3047       INTEGER(iwp)        ::  kc                    !<
3048
3049       REAL(wp), PARAMETER ::  cfw = 0.2_wp          !<
3050       REAL(wp), PARAMETER ::  c_tkef = 0.6_wp       !<
3051       REAL(wp)            ::  fw                    !<
3052       REAL(wp), PARAMETER ::  fw0 = 0.9_wp          !<
3053       REAL(wp)            ::  glsf                  !<
3054       REAL(wp)            ::  glsc                  !<
3055       REAL(wp)            ::  height                !<
3056       REAL(wp), PARAMETER ::  p13 = 1.0_wp/3.0_wp   !<
3057       REAL(wp), PARAMETER ::  p23 = 2.0_wp/3.0_wp   !<       
3058
3059!
3060       IF ( .NOT. rans_mode  .AND.  .NOT. rans_mode_parent )  THEN
3061          IF ( nest_bound_l )  THEN
3062             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
3063             tkefactor_l = 0.0_wp
3064             i = nxl - 1
3065             DO  j = nysg, nyng
3066                k_wall = get_topography_top_index_ji( j, i, 's' )
3067
3068                DO  k = k_wall + 1, nzt
3069                   kc     = kco(k) + 1
3070                   glsf   = ( dx * dy * dzu(k) )**p13
3071                   glsc   = ( cg%dx * cg%dy *cg%dzu(kc) )**p13
3072                   height = zu(k) - zu(k_wall)
3073                   fw     = EXP( -cfw * height / glsf )
3074                   tkefactor_l(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
3075                                                 ( glsf / glsc )**p23 )
3076                ENDDO
3077                tkefactor_l(k_wall,j) = c_tkef * fw0
3078             ENDDO
3079          ENDIF
3080
3081          IF ( nest_bound_r )  THEN
3082             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
3083             tkefactor_r = 0.0_wp
3084             i = nxr + 1
3085             DO  j = nysg, nyng
3086                k_wall = get_topography_top_index_ji( j, i, 's' )
3087
3088                DO  k = k_wall + 1, nzt
3089                   kc     = kco(k) + 1
3090                   glsf   = ( dx * dy * dzu(k) )**p13
3091                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
3092                   height = zu(k) - zu(k_wall)
3093                   fw     = EXP( -cfw * height / glsf )
3094                   tkefactor_r(k,j) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
3095                                                 ( glsf / glsc )**p23 )
3096                ENDDO
3097                tkefactor_r(k_wall,j) = c_tkef * fw0
3098             ENDDO
3099          ENDIF
3100
3101          IF ( nest_bound_s )  THEN
3102             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
3103             tkefactor_s = 0.0_wp
3104             j = nys - 1
3105             DO  i = nxlg, nxrg
3106                k_wall = get_topography_top_index_ji( j, i, 's' )
3107               
3108                DO  k = k_wall + 1, nzt
3109   
3110                   kc     = kco(k) + 1
3111                   glsf   = ( dx * dy * dzu(k) )**p13
3112                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) ) ** p13
3113                   height = zu(k) - zu(k_wall)
3114                   fw     = EXP( -cfw*height / glsf )
3115                   tkefactor_s(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *      &
3116                        ( glsf / glsc )**p23 )
3117                ENDDO
3118                tkefactor_s(k_wall,i) = c_tkef * fw0
3119             ENDDO
3120          ENDIF
3121
3122          IF ( nest_bound_n )  THEN
3123             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
3124             tkefactor_n = 0.0_wp
3125             j = nyn + 1
3126             DO  i = nxlg, nxrg
3127                k_wall = get_topography_top_index_ji( j, i, 's' )
3128
3129                DO  k = k_wall + 1, nzt
3130
3131                   kc     = kco(k) + 1
3132                   glsf   = ( dx * dy * dzu(k) )**p13
3133                   glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
3134                   height = zu(k) - zu(k_wall)
3135                   fw     = EXP( -cfw * height / glsf )
3136                   tkefactor_n(k,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *     &
3137                                                 ( glsf / glsc )**p23 )
3138                ENDDO
3139                tkefactor_n(k_wall,i) = c_tkef * fw0
3140             ENDDO
3141          ENDIF
3142
3143          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
3144          k = nzt
3145
3146          DO  i = nxlg, nxrg
3147             DO  j = nysg, nyng
3148!
3149!--             Determine vertical index for local topography top
3150                k_wall = get_topography_top_index_ji( j, i, 's' )
3151
3152                kc     = kco(k) + 1
3153                glsf   = ( dx * dy * dzu(k) )**p13
3154                glsc   = ( cg%dx * cg%dy * cg%dzu(kc) )**p13
3155                height = zu(k) - zu(k_wall)
3156                fw     = EXP( -cfw * height / glsf )
3157                tkefactor_t(j,i) = c_tkef * ( fw0 * fw + ( 1.0_wp - fw ) *        &
3158                                              ( glsf / glsc )**p23 )
3159             ENDDO
3160          ENDDO
3161!
3162!--    RANS mode
3163       ELSE
3164          IF ( nest_bound_l )  THEN
3165             ALLOCATE( tkefactor_l(nzb:nzt+1,nysg:nyng) )
3166             tkefactor_l = 1.0_wp
3167          ENDIF
3168          IF ( nest_bound_r )  THEN
3169             ALLOCATE( tkefactor_r(nzb:nzt+1,nysg:nyng) )
3170             tkefactor_r = 1.0_wp
3171          ENDIF
3172          IF ( nest_bound_s )  THEN
3173             ALLOCATE( tkefactor_s(nzb:nzt+1,nxlg:nxrg) )
3174             tkefactor_s = 1.0_wp
3175          ENDIF
3176          IF ( nest_bound_n )  THEN
3177             ALLOCATE( tkefactor_n(nzb:nzt+1,nxlg:nxrg) )
3178             tkefactor_n = 1.0_wp
3179          ENDIF
3180
3181          ALLOCATE( tkefactor_t(nysg:nyng,nxlg:nxrg) )
3182          tkefactor_t = 1.0_wp
3183
3184       ENDIF
3185     
3186    END SUBROUTINE pmci_init_tkefactor
3187
3188#endif
3189 END SUBROUTINE pmci_setup_child
3190
3191
3192
3193 SUBROUTINE pmci_setup_coordinates
3194
3195#if defined( __parallel )
3196    IMPLICIT NONE
3197
3198    INTEGER(iwp) ::  i   !<
3199    INTEGER(iwp) ::  j   !<
3200
3201!
3202!-- Create coordinate arrays.
3203    ALLOCATE( coord_x(-nbgp:nx+nbgp) )
3204    ALLOCATE( coord_y(-nbgp:ny+nbgp) )
3205     
3206    DO  i = -nbgp, nx + nbgp
3207       coord_x(i) = lower_left_coord_x + i * dx
3208    ENDDO
3209     
3210    DO  j = -nbgp, ny + nbgp
3211       coord_y(j) = lower_left_coord_y + j * dy
3212    ENDDO
3213
3214#endif
3215 END SUBROUTINE pmci_setup_coordinates
3216
3217
3218
3219
3220 SUBROUTINE pmci_set_array_pointer( name, child_id, nz_cl, n )
3221
3222    IMPLICIT NONE
3223
3224    INTEGER(iwp), INTENT(IN)          ::  child_id    !<
3225    INTEGER(iwp), INTENT(IN)          ::  nz_cl       !<
3226    INTEGER(iwp), INTENT(IN),OPTIONAL ::  n           !< index of chemical species
3227
3228    CHARACTER(LEN=*), INTENT(IN) ::  name        !<
3229
3230#if defined( __parallel )
3231    INTEGER(iwp) ::  ierr        !<
3232    INTEGER(iwp) ::  istat       !<
3233
3234    REAL(wp), POINTER, DIMENSION(:,:)     ::  p_2d        !<
3235    REAL(wp), POINTER, DIMENSION(:,:)     ::  p_2d_sec    !<
3236    REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d        !<
3237    REAL(wp), POINTER, DIMENSION(:,:,:)   ::  p_3d_sec    !<
3238    INTEGER(idp), POINTER, DIMENSION(:,:) ::  i_2d        !<
3239
3240
3241    NULLIFY( p_3d )
3242    NULLIFY( p_2d )
3243    NULLIFY( i_2d )
3244
3245!
3246!-- List of array names, which can be coupled.
3247!-- In case of 3D please change also the second array for the pointer version
3248    IF ( TRIM(name) == "u"          )  p_3d => u
3249    IF ( TRIM(name) == "v"          )  p_3d => v
3250    IF ( TRIM(name) == "w"          )  p_3d => w
3251    IF ( TRIM(name) == "e"          )  p_3d => e
3252    IF ( TRIM(name) == "pt"         )  p_3d => pt
3253    IF ( TRIM(name) == "q"          )  p_3d => q
3254    IF ( TRIM(name) == "qc"         )  p_3d => qc
3255    IF ( TRIM(name) == "qr"         )  p_3d => qr
3256    IF ( TRIM(name) == "nr"         )  p_3d => nr
3257    IF ( TRIM(name) == "nc"         )  p_3d => nc
3258    IF ( TRIM(name) == "s"          )  p_3d => s
3259    IF ( TRIM(name) == "diss"       )  p_3d => diss   
3260    IF ( TRIM(name) == "nr_part"    )   i_2d => nr_part
3261    IF ( TRIM(name) == "part_adr"   )  i_2d => part_adr
3262    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d => chem_species(n)%conc
3263
3264!
3265!-- Next line is just an example for a 2D array (not active for coupling!)
3266!-- Please note, that z0 has to be declared as TARGET array in modules.f90
3267!    IF ( TRIM(name) == "z0" )    p_2d => z0
3268
3269#if defined( __nopointer )
3270    IF ( ASSOCIATED( p_3d ) )  THEN
3271       CALL pmc_s_set_dataarray( child_id, p_3d, nz_cl, nz )
3272    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3273       CALL pmc_s_set_dataarray( child_id, p_2d )
3274    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
3275       CALL pmc_s_set_dataarray( child_id, i_2d )
3276    ELSE
3277!
3278!--    Give only one message for the root domain
3279       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
3280
3281          message_string = 'pointer for array "' // TRIM( name ) //            &
3282                           '" can''t be associated'
3283          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
3284       ELSE
3285!
3286!--       Avoid others to continue
3287          CALL MPI_BARRIER( comm2d, ierr )
3288       ENDIF
3289    ENDIF
3290#else
3291    IF ( TRIM(name) == "u"    )  p_3d_sec => u_2
3292    IF ( TRIM(name) == "v"    )  p_3d_sec => v_2
3293    IF ( TRIM(name) == "w"    )  p_3d_sec => w_2
3294    IF ( TRIM(name) == "e"    )  p_3d_sec => e_2
3295    IF ( TRIM(name) == "pt"   )  p_3d_sec => pt_2
3296    IF ( TRIM(name) == "q"    )  p_3d_sec => q_2
3297    IF ( TRIM(name) == "qc"   )  p_3d_sec => qc_2
3298    IF ( TRIM(name) == "qr"   )  p_3d_sec => qr_2
3299    IF ( TRIM(name) == "nr"   )  p_3d_sec => nr_2
3300    IF ( TRIM(name) == "nc"   )  p_3d_sec => nc_2
3301    IF ( TRIM(name) == "s"    )  p_3d_sec => s_2
3302    IF ( TRIM(name) == "diss" )  p_3d_sec => diss_2
3303    IF ( INDEX( TRIM(name), "chem_" ) /= 0 )  p_3d_sec => spec_conc_2(:,:,:,n)
3304
3305    IF ( ASSOCIATED( p_3d ) )  THEN
3306       CALL pmc_s_set_dataarray( child_id, p_3d, nz_cl, nz,                    &
3307                                 array_2 = p_3d_sec )
3308    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3309       CALL pmc_s_set_dataarray( child_id, p_2d )
3310    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
3311       CALL pmc_s_set_dataarray( child_id, i_2d )
3312    ELSE
3313!
3314!--    Give only one message for the root domain
3315       IF ( myid == 0  .AND.  cpl_id == 1 )  THEN
3316
3317          message_string = 'pointer for array "' // TRIM( name ) //            &
3318                           '" can''t be associated'
3319          CALL message( 'pmci_set_array_pointer', 'PA0117', 3, 2, 0, 6, 0 )
3320       ELSE
3321!
3322!--       Avoid others to continue
3323          CALL MPI_BARRIER( comm2d, ierr )
3324       ENDIF
3325
3326   ENDIF
3327#endif
3328
3329#endif
3330END SUBROUTINE pmci_set_array_pointer
3331
3332INTEGER FUNCTION get_number_of_childs ()
3333   IMPLICIT NONE
3334
3335   get_number_of_childs = SIZE( pmc_parent_for_child ) - 1
3336
3337   RETURN
3338END FUNCTION get_number_of_childs
3339
3340INTEGER FUNCTION get_childid (id_index)
3341   IMPLICIT NONE
3342
3343   INTEGER,INTENT(IN)                 :: id_index
3344
3345   get_childid = pmc_parent_for_child(id_index)
3346
3347   RETURN
3348END FUNCTION get_childid
3349
3350SUBROUTINE  get_child_edges (m, lx_coord, lx_coord_b, rx_coord, rx_coord_b,    &
3351                               sy_coord, sy_coord_b, ny_coord, ny_coord_b,     &
3352                               uz_coord, uz_coord_b)
3353   IMPLICIT NONE
3354   INTEGER,INTENT(IN)             ::  m
3355   REAL(wp),INTENT(OUT)           ::  lx_coord, lx_coord_b
3356   REAL(wp),INTENT(OUT)           ::  rx_coord, rx_coord_b
3357   REAL(wp),INTENT(OUT)           ::  sy_coord, sy_coord_b
3358   REAL(wp),INTENT(OUT)           ::  ny_coord, ny_coord_b
3359   REAL(wp),INTENT(OUT)           ::  uz_coord, uz_coord_b
3360
3361   lx_coord = childgrid(m)%lx_coord
3362   rx_coord = childgrid(m)%rx_coord
3363   sy_coord = childgrid(m)%sy_coord
3364   ny_coord = childgrid(m)%ny_coord
3365   uz_coord = childgrid(m)%uz_coord
3366
3367   lx_coord_b = childgrid(m)%lx_coord_b
3368   rx_coord_b = childgrid(m)%rx_coord_b
3369   sy_coord_b = childgrid(m)%sy_coord_b
3370   ny_coord_b = childgrid(m)%ny_coord_b
3371   uz_coord_b = childgrid(m)%uz_coord_b
3372
3373END SUBROUTINE get_child_edges
3374
3375SUBROUTINE  get_child_gridspacing (m, dx,dy,dz)
3376
3377   IMPLICIT NONE
3378   INTEGER,INTENT(IN)             ::  m
3379   REAL(wp),INTENT(OUT)           ::  dx,dy
3380   REAL(wp),INTENT(OUT),OPTIONAL  ::  dz
3381
3382   dx = childgrid(m)%dx
3383   dy = childgrid(m)%dy
3384   IF(PRESENT(dz))   THEN
3385      dz = childgrid(m)%dz
3386   ENDIF
3387
3388END SUBROUTINE get_child_gridspacing
3389
3390SUBROUTINE pmci_create_child_arrays( name, is, ie, js, je, nzc,)
3391
3392    IMPLICIT NONE
3393
3394    CHARACTER(LEN=*), INTENT(IN) ::  name    !<
3395
3396    INTEGER(iwp), INTENT(IN) ::  ie      !<
3397    INTEGER(iwp), INTENT(IN) ::  is      !<
3398    INTEGER(iwp), INTENT(IN) ::  je      !<
3399    INTEGER(iwp), INTENT(IN) ::  js      !<
3400    INTEGER(iwp), INTENT(IN) ::  nzc     !<  Note that nzc is cg%nz
3401
3402    INTEGER(iwp), INTENT(IN), OPTIONAL ::  n  !< number of chemical species
3403
3404#if defined( __parallel )
3405    INTEGER(iwp) ::  ierr    !<
3406    INTEGER(iwp) ::  istat   !<
3407
3408    REAL(wp), POINTER,DIMENSION(:,:)       ::  p_2d    !<
3409    REAL(wp), POINTER,DIMENSION(:,:,:)     ::  p_3d    !<
3410    INTEGER(idp), POINTER,DIMENSION(:,:)   ::  i_2d    !<
3411
3412
3413    NULLIFY( p_3d )
3414    NULLIFY( p_2d )
3415    NULLIFY( i_2d )
3416
3417!
3418!-- List of array names, which can be coupled
3419    IF ( TRIM( name ) == "u" )  THEN
3420       IF ( .NOT. ALLOCATED( uc ) )  ALLOCATE( uc(0:nzc+1,js:je,is:ie) )
3421       p_3d => uc
3422    ELSEIF ( TRIM( name ) == "v" )  THEN
3423       IF ( .NOT. ALLOCATED( vc ) )  ALLOCATE( vc(0:nzc+1,js:je,is:ie) )
3424       p_3d => vc
3425    ELSEIF ( TRIM( name ) == "w" )  THEN
3426       IF ( .NOT. ALLOCATED( wc ) )  ALLOCATE( wc(0:nzc+1,js:je,is:ie) )
3427       p_3d => wc
3428    ELSEIF ( TRIM( name ) == "e" )  THEN
3429       IF ( .NOT. ALLOCATED( ec ) )  ALLOCATE( ec(0:nzc+1,js:je,is:ie) )
3430       p_3d => ec
3431    ELSEIF ( TRIM( name ) == "diss" )  THEN
3432       IF ( .NOT. ALLOCATED( dissc ) )  ALLOCATE( dissc(0:nzc+1,js:je,is:ie) )
3433       p_3d => dissc
3434    ELSEIF ( TRIM( name ) == "pt")  THEN
3435       IF ( .NOT. ALLOCATED( ptc ) )  ALLOCATE( ptc(0:nzc+1,js:je,is:ie) )
3436       p_3d => ptc
3437    ELSEIF ( TRIM( name ) == "q")  THEN
3438       IF ( .NOT. ALLOCATED( q_c ) )  ALLOCATE( q_c(0:nzc+1,js:je,is:ie) )
3439       p_3d => q_c
3440    ELSEIF ( TRIM( name ) == "qc")  THEN
3441       IF ( .NOT. ALLOCATED( qcc ) )  ALLOCATE( qcc(0:nzc+1,js:je,is:ie) )
3442       p_3d => qcc
3443    ELSEIF ( TRIM( name ) == "qr")  THEN
3444       IF ( .NOT. ALLOCATED( qrc ) )  ALLOCATE( qrc(0:nzc+1,js:je,is:ie) )
3445       p_3d => qrc
3446    ELSEIF ( TRIM( name ) == "nr")  THEN
3447       IF ( .NOT. ALLOCATED( nrc ) )  ALLOCATE( nrc(0:nzc+1,js:je,is:ie) )
3448       p_3d => nrc
3449    ELSEIF ( TRIM( name ) == "nc")  THEN
3450       IF ( .NOT. ALLOCATED( ncc ) )  ALLOCATE( ncc(0:nzc+1,js:je,is:ie) )
3451       p_3d => ncc
3452    ELSEIF ( TRIM( name ) == "s")  THEN
3453       IF ( .NOT. ALLOCATED( sc ) )  ALLOCATE( sc(0:nzc+1,js:je,is:ie) )
3454       p_3d => sc
3455    ELSEIF ( TRIM( name ) == "nr_part") THEN
3456       IF ( .NOT. ALLOCATED( nr_partc ) )  ALLOCATE( nr_partc(js:je,is:ie) )
3457       i_2d => nr_partc
3458    ELSEIF ( TRIM( name ) == "part_adr") THEN
3459       IF ( .NOT. ALLOCATED( part_adrc ) )  ALLOCATE( part_adrc(js:je,is:ie) )
3460       i_2d => part_adrc
3461    ELSEIF ( TRIM( name(1:5) ) == "chem_" )  THEN
3462       IF ( .NOT. ALLOCATED( chem_spec_c ) )                                   &
3463          ALLOCATE( chem_spec_c(0:nzc+1,js:je,is:ie,1:nspec) )
3464       p_3d => chem_spec_c(:,:,:,n)
3465    !ELSEIF (trim(name) == "z0") then
3466       !IF (.not.allocated(z0c))  allocate(z0c(js:je, is:ie))
3467       !p_2d => z0c
3468    ENDIF
3469
3470    IF ( ASSOCIATED( p_3d ) )  THEN
3471       CALL pmc_c_set_dataarray( p_3d )
3472    ELSEIF ( ASSOCIATED( p_2d ) )  THEN
3473       CALL pmc_c_set_dataarray( p_2d )
3474    ELSEIF ( ASSOCIATED( i_2d ) )  THEN
3475       CALL pmc_c_set_dataarray( i_2d )
3476    ELSE
3477!
3478!--    Give only one message for the first child domain
3479       IF ( myid == 0  .AND.  cpl_id == 2 )  THEN
3480
3481          message_string = 'pointer for array "' // TRIM( name ) //            &
3482                           '" can''t be associated'
3483          CALL message( 'pmci_create_child_arrays', 'PA0170', 3, 2, 0, 6, 0 )
3484       ELSE
3485!
3486!--       Prevent others from continuing
3487          CALL MPI_BARRIER( comm2d, ierr )
3488       ENDIF
3489    ENDIF
3490
3491#endif
3492 END SUBROUTINE pmci_create_child_arrays
3493
3494
3495
3496 SUBROUTINE pmci_parent_initialize
3497
3498!
3499!-- Send data for the children in order to let them create initial
3500!-- conditions by interpolating the parent-domain fields.
3501#if defined( __parallel )
3502    IMPLICIT NONE
3503
3504    INTEGER(iwp) ::  child_id    !<
3505    INTEGER(iwp) ::  m           !<
3506
3507    REAL(wp) ::  waittime        !<
3508
3509
3510    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
3511       child_id = pmc_parent_for_child(m)
3512       CALL pmc_s_fillbuffer( child_id, waittime=waittime )
3513    ENDDO
3514
3515#endif
3516 END SUBROUTINE pmci_parent_initialize
3517
3518
3519
3520 SUBROUTINE pmci_child_initialize
3521
3522!
3523!-- Create initial conditions for the current child domain by interpolating
3524!-- the parent-domain fields.
3525#if defined( __parallel )
3526    IMPLICIT NONE
3527
3528    INTEGER(iwp) ::  i          !<
3529    INTEGER(iwp) ::  icl        !<
3530    INTEGER(iwp) ::  icr        !<
3531    INTEGER(iwp) ::  j          !<
3532    INTEGER(iwp) ::  jcn        !<
3533    INTEGER(iwp) ::  jcs        !<
3534    INTEGER(iwp) ::  k          !<
3535    INTEGER(iwp) ::  n          !< running index for chemical species
3536
3537    REAL(wp) ::  waittime       !<
3538
3539!
3540!-- Root model is never anyone's child
3541    IF ( cpl_id > 1 )  THEN
3542!
3543!--    Child domain boundaries in the parent index space
3544       icl = coarse_bound(1)
3545       icr = coarse_bound(2)
3546       jcs = coarse_bound(3)
3547       jcn = coarse_bound(4)
3548!
3549!--    Get data from the parent
3550       CALL pmc_c_getbuffer( waittime = waittime )
3551!
3552!--    The interpolation.
3553       CALL pmci_interp_tril_all ( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,   &
3554                                   r2yo, r1zo, r2zo, 'u' )
3555       CALL pmci_interp_tril_all ( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,   &
3556                                   r2yv, r1zo, r2zo, 'v' )
3557       CALL pmci_interp_tril_all ( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,   &
3558                                   r2yo, r1zw, r2zw, 'w' )
3559
3560       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.          &
3561            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.           &
3562               .NOT. constant_diffusion ) )  THEN
3563          CALL pmci_interp_tril_all ( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo, &
3564                                      r2yo, r1zo, r2zo, 'e' )
3565       ENDIF
3566
3567       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
3568          CALL pmci_interp_tril_all ( diss,  dissc,  ico, jco, kco, r1xo, r2xo,&
3569                                      r1yo, r2yo, r1zo, r2zo, 's' )
3570       ENDIF
3571
3572       IF ( .NOT. neutral )  THEN
3573          CALL pmci_interp_tril_all ( pt, ptc, ico, jco, kco, r1xo, r2xo,      &
3574                                      r1yo, r2yo, r1zo, r2zo, 's' )
3575       ENDIF
3576
3577       IF ( humidity )  THEN
3578
3579          CALL pmci_interp_tril_all ( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo, &
3580                                      r2yo, r1zo, r2zo, 's' )
3581
3582          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
3583             CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,   &
3584                                          r1yo, r2yo, r1zo, r2zo, 's' ) 
3585             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,   &
3586                                         r1yo, r2yo, r1zo, r2zo, 's' )   
3587          ENDIF
3588
3589          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
3590             CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,   &
3591                                         r1yo, r2yo, r1zo, r2zo, 's' )
3592             CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,   &
3593                                         r1yo, r2yo, r1zo, r2zo, 's' )
3594          ENDIF
3595
3596       ENDIF
3597
3598       IF ( passive_scalar )  THEN
3599          CALL pmci_interp_tril_all ( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,  &
3600                                      r2yo, r1zo, r2zo, 's' )
3601       ENDIF
3602
3603       IF ( air_chemistry )  THEN
3604          DO  n = 1, nspec
3605             CALL pmci_interp_tril_all ( chem_species(n)%conc,                 &
3606                                         chem_spec_c(:,:,:,n),                 &
3607                                         ico, jco, kco, r1xo, r2xo, r1yo,      &
3608                                         r2yo, r1zo, r2zo, 's' )
3609          ENDDO
3610       ENDIF
3611
3612       IF ( topography /= 'flat' )  THEN
3613!
3614!--       Inside buildings set velocities and TKE back to zero.
3615!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
3616!--       maybe revise later.
3617          DO   i = nxlg, nxrg
3618             DO   j = nysg, nyng
3619                DO  k = nzb, nzt
3620                   u(k,j,i)   = MERGE( u(k,j,i), 0.0_wp,                       &
3621                                       BTEST( wall_flags_0(k,j,i), 1 ) )
3622                   v(k,j,i)   = MERGE( v(k,j,i), 0.0_wp,                       &
3623                                       BTEST( wall_flags_0(k,j,i), 2 ) )
3624                   w(k,j,i)   = MERGE( w(k,j,i), 0.0_wp,                       &
3625                                       BTEST( wall_flags_0(k,j,i), 3 ) )
3626!                    e(k,j,i)   = MERGE( e(k,j,i), 0.0_wp,                       &
3627!                                        BTEST( wall_flags_0(k,j,i), 0 ) )
3628                   u_p(k,j,i) = MERGE( u_p(k,j,i), 0.0_wp,                     &
3629                                       BTEST( wall_flags_0(k,j,i), 1 ) )
3630                   v_p(k,j,i) = MERGE( v_p(k,j,i), 0.0_wp,                     &
3631                                       BTEST( wall_flags_0(k,j,i), 2 ) )
3632                   w_p(k,j,i) = MERGE( w_p(k,j,i), 0.0_wp,                     &
3633                                       BTEST( wall_flags_0(k,j,i), 3 ) )
3634!                    e_p(k,j,i) = MERGE( e_p(k,j,i), 0.0_wp,                     &
3635!                                        BTEST( wall_flags_0(k,j,i), 0 ) )
3636                ENDDO
3637             ENDDO
3638          ENDDO
3639       ENDIF
3640    ENDIF
3641
3642
3643 CONTAINS
3644
3645
3646    SUBROUTINE pmci_interp_tril_all( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y,    &
3647                                     r1z, r2z, var )
3648!
3649!--    Interpolation of the internal values for the child-domain initialization
3650!--    This subroutine is based on trilinear interpolation.
3651!--    Coding based on interp_tril_lr/sn/t
3652       IMPLICIT NONE
3653
3654       CHARACTER(LEN=1), INTENT(IN) :: var  !<
3655
3656       INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic    !<
3657       INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc    !<
3658       INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc    !<
3659
3660       INTEGER(iwp) ::  i        !<
3661       INTEGER(iwp) ::  ib       !<
3662       INTEGER(iwp) ::  ie       !<
3663       INTEGER(iwp) ::  j        !<
3664       INTEGER(iwp) ::  jb       !<
3665       INTEGER(iwp) ::  je       !<
3666       INTEGER(iwp) ::  k        !<
3667       INTEGER(iwp) ::  k_wall   !<
3668       INTEGER(iwp) ::  k1       !<
3669       INTEGER(iwp) ::  kb       !<
3670       INTEGER(iwp) ::  kbc      !<
3671       INTEGER(iwp) ::  l        !<
3672       INTEGER(iwp) ::  m        !<
3673       INTEGER(iwp) ::  n        !<
3674
3675       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(INOUT) :: f !<
3676       REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr), INTENT(IN) :: fc       !<
3677       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r1x   !<
3678       REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN) :: r2x   !<
3679       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r1y   !<
3680       REAL(wp), DIMENSION(nysg:nyng), INTENT(IN) :: r2y   !<
3681       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r1z   !<
3682       REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN) :: r2z   !<
3683
3684       REAL(wp) ::  fk         !<
3685       REAL(wp) ::  fkj        !<
3686       REAL(wp) ::  fkjp       !<
3687       REAL(wp) ::  fkp        !<
3688       REAL(wp) ::  fkpj       !<
3689       REAL(wp) ::  fkpjp      !<
3690       REAL(wp) ::  logratio   !<
3691       REAL(wp) ::  logzuc1    !<
3692       REAL(wp) ::  zuc1       !<
3693       REAL(wp) ::  z0_topo    !<  roughness at vertical walls
3694
3695
3696       ib = nxl
3697       ie = nxr
3698       jb = nys
3699       je = nyn
3700       IF ( nesting_mode /= 'vertical' )  THEN
3701          IF ( nest_bound_l )  THEN
3702             ib = nxl - 1
3703!
3704!--          For u, nxl is a ghost node, but not for the other variables
3705             IF ( var == 'u' )  THEN
3706                ib = nxl
3707             ENDIF
3708          ENDIF
3709          IF ( nest_bound_s )  THEN
3710             jb = nys - 1
3711!
3712!--          For v, nys is a ghost node, but not for the other variables
3713             IF ( var == 'v' )  THEN
3714                jb = nys
3715             ENDIF
3716          ENDIF
3717          IF ( nest_bound_r )  THEN
3718             ie = nxr + 1
3719          ENDIF
3720          IF ( nest_bound_n )  THEN
3721             je = nyn + 1
3722          ENDIF
3723       ENDIF
3724!
3725!--    Trilinear interpolation.
3726       DO  i = ib, ie
3727          DO  j = jb, je
3728!
3729!--          Determine the vertical index of the first node above the
3730!--          topography top at grid point (j,i) in order to not overwrite
3731!--          the bottom BC.
3732             kb = get_topography_top_index_ji( j, i, TRIM ( var ) ) + 1
3733             DO  k = kb, nzt + 1
3734                l = ic(i)
3735                m = jc(j)
3736                n = kc(k)
3737                fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
3738                fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
3739                fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
3740                fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
3741                fk       = r1y(j) * fkj  + r2y(j) * fkjp
3742                fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
3743                f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
3744             ENDDO
3745          ENDDO
3746       ENDDO
3747!
3748!--    Correct the interpolated values of u and v in near-wall nodes, i.e. in
3749!--    the nodes below the coarse-grid nodes with k=1. The corrction is only
3750!--    made over horizontal wall surfaces in this phase. For the nest boundary
3751!--    conditions, a corresponding correction is made for all vertical walls,
3752!--    too.
3753       IF ( constant_flux_layer .AND. ( var == 'u' .OR. var == 'v' ) )  THEN
3754          z0_topo = roughness_length
3755          DO  i = ib, nxr
3756             DO  j = jb, nyn
3757!
3758!--             Determine vertical index of topography top at grid point (j,i)
3759                k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
3760!
3761!--             kbc is the first coarse-grid point above the surface
3762                kbc = 1
3763                DO  WHILE ( cg%zu(kbc) < zu(k_wall) )
3764                   kbc = kbc + 1
3765                ENDDO
3766                zuc1 = cg%zu(kbc)
3767                k1   = k_wall + 1
3768                DO  WHILE ( zu(k1) < zuc1 )
3769                   k1 = k1 + 1
3770                ENDDO
3771                logzuc1 = LOG( ( zu(k1) - zu(k_wall) ) / z0_topo )
3772
3773                k = k_wall + 1
3774                DO  WHILE ( zu(k) < zuc1 )
3775                   logratio = ( LOG( ( zu(k) - zu(k_wall) ) / z0_topo ) ) /    &
3776                                logzuc1
3777                   f(k,j,i) = logratio * f(k1,j,i)
3778                   k  = k + 1
3779                ENDDO
3780                f(k_wall,j,i) = 0.0_wp
3781             ENDDO
3782          ENDDO
3783
3784       ELSEIF ( var == 'w' )  THEN
3785
3786          DO  i = ib, nxr
3787              DO  j = jb, nyn
3788!
3789!--              Determine vertical index of topography top at grid point (j,i)
3790                 k_wall = get_topography_top_index_ji( j, i, 'w' )
3791
3792                 f(k_wall,j,i) = 0.0_wp
3793              ENDDO
3794           ENDDO
3795
3796       ENDIF
3797
3798    END SUBROUTINE pmci_interp_tril_all
3799
3800#endif
3801 END SUBROUTINE pmci_child_initialize
3802
3803
3804
3805 SUBROUTINE pmci_check_setting_mismatches
3806!
3807!-- Check for mismatches between settings of master and child variables
3808!-- (e.g., all children have to follow the end_time settings of the root model).
3809!-- The root model overwrites variables in the other models, so these variables
3810!-- only need to be set once in file PARIN.
3811
3812#if defined( __parallel )
3813
3814    USE control_parameters,                                                    &
3815        ONLY:  dt_restart, end_time, message_string, restart_time, time_restart
3816
3817    IMPLICIT NONE
3818
3819    INTEGER ::  ierr
3820
3821    REAL(wp) ::  dt_restart_root
3822    REAL(wp) ::  end_time_root
3823    REAL(wp) ::  restart_time_root
3824    REAL(wp) ::  time_restart_root
3825
3826!
3827!-- Check the time to be simulated.
3828!-- Here, and in the following, the root process communicates the respective
3829!-- variable to all others, and its value will then be compared with the local
3830!-- values.
3831    IF ( pmc_is_rootmodel() )  end_time_root = end_time
3832    CALL MPI_BCAST( end_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3833
3834    IF ( .NOT. pmc_is_rootmodel() )  THEN
3835       IF ( end_time /= end_time_root )  THEN
3836          WRITE( message_string, * )  'mismatch between root model and ',      &
3837               'child settings &   end_time(root) = ', end_time_root,          &
3838               ' &   end_time(child) = ', end_time, ' & child value is set',   &
3839               ' to root value'
3840          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3841                        0 )
3842          end_time = end_time_root
3843       ENDIF
3844    ENDIF
3845!
3846!-- Same for restart time
3847    IF ( pmc_is_rootmodel() )  restart_time_root = restart_time
3848    CALL MPI_BCAST( restart_time_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3849
3850    IF ( .NOT. pmc_is_rootmodel() )  THEN
3851       IF ( restart_time /= restart_time_root )  THEN
3852          WRITE( message_string, * )  'mismatch between root model and ',      &
3853               'child settings &   restart_time(root) = ', restart_time_root,  &
3854               ' &   restart_time(child) = ', restart_time, ' & child ',       &
3855               'value is set to root value'
3856          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3857                        0 )
3858          restart_time = restart_time_root
3859       ENDIF
3860    ENDIF
3861!
3862!-- Same for dt_restart
3863    IF ( pmc_is_rootmodel() )  dt_restart_root = dt_restart
3864    CALL MPI_BCAST( dt_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3865
3866    IF ( .NOT. pmc_is_rootmodel() )  THEN
3867       IF ( dt_restart /= dt_restart_root )  THEN
3868          WRITE( message_string, * )  'mismatch between root model and ',      &
3869               'child settings &   dt_restart(root) = ', dt_restart_root,      &
3870               ' &   dt_restart(child) = ', dt_restart, ' & child ',           &
3871               'value is set to root value'
3872          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3873                        0 )
3874          dt_restart = dt_restart_root
3875       ENDIF
3876    ENDIF
3877!
3878!-- Same for time_restart
3879    IF ( pmc_is_rootmodel() )  time_restart_root = time_restart
3880    CALL MPI_BCAST( time_restart_root, 1, MPI_REAL, 0, comm_world_nesting, ierr )
3881
3882    IF ( .NOT. pmc_is_rootmodel() )  THEN
3883       IF ( time_restart /= time_restart_root )  THEN
3884          WRITE( message_string, * )  'mismatch between root model and ',      &
3885               'child settings &   time_restart(root) = ', time_restart_root,  &
3886               ' &   time_restart(child) = ', time_restart, ' & child ',       &
3887               'value is set to root value'
3888          CALL message( 'pmci_check_setting_mismatches', 'PA0419', 0, 1, 0, 6, &
3889                        0 )
3890          time_restart = time_restart_root
3891       ENDIF
3892    ENDIF
3893
3894#endif
3895
3896 END SUBROUTINE pmci_check_setting_mismatches
3897
3898
3899
3900 SUBROUTINE pmci_ensure_nest_mass_conservation
3901
3902!
3903!-- Adjust the volume-flow rate through the top boundary so that the net volume
3904!-- flow through all boundaries of the current nest domain becomes zero.
3905    IMPLICIT NONE
3906
3907    INTEGER(iwp) ::  i                           !<
3908    INTEGER(iwp) ::  ierr                        !<
3909    INTEGER(iwp) ::  j                           !<
3910    INTEGER(iwp) ::  k                           !<
3911
3912    REAL(wp) ::  dxdy                            !<
3913    REAL(wp) ::  innor                           !<
3914    REAL(wp) ::  w_lt                            !<
3915    REAL(wp), DIMENSION(1:3) ::  volume_flow_l   !<
3916
3917!
3918!-- Sum up the volume flow through the left/right boundaries
3919    volume_flow(1)   = 0.0_wp
3920    volume_flow_l(1) = 0.0_wp
3921
3922    IF ( nest_bound_l )  THEN
3923       i = 0
3924       innor = dy
3925       DO   j = nys, nyn
3926          DO   k = nzb+1, nzt
3927             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
3928                                     * MERGE( 1.0_wp, 0.0_wp,                  &
3929                                              BTEST( wall_flags_0(k,j,i), 1 ) )
3930          ENDDO
3931       ENDDO
3932    ENDIF
3933
3934    IF ( nest_bound_r )  THEN
3935       i = nx + 1
3936       innor = -dy
3937       DO   j = nys, nyn
3938          DO   k = nzb+1, nzt
3939             volume_flow_l(1) = volume_flow_l(1) + innor * u(k,j,i) * dzw(k)   &
3940                                     * MERGE( 1.0_wp, 0.0_wp,                  &
3941                                              BTEST( wall_flags_0(k,j,i), 1 ) )
3942          ENDDO
3943       ENDDO
3944    ENDIF
3945
3946#if defined( __parallel )
3947    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3948    CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL,         &
3949                        MPI_SUM, comm2d, ierr )
3950#else
3951    volume_flow(1) = volume_flow_l(1)
3952#endif
3953   
3954!
3955!-- Sum up the volume flow through the south/north boundaries
3956    volume_flow(2)   = 0.0_wp
3957    volume_flow_l(2) = 0.0_wp
3958
3959    IF ( nest_bound_s )  THEN
3960       j = 0
3961       innor = dx
3962       DO   i = nxl, nxr
3963          DO   k = nzb+1, nzt
3964             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
3965                                     * MERGE( 1.0_wp, 0.0_wp,                  &
3966                                              BTEST( wall_flags_0(k,j,i), 2 ) )
3967          ENDDO
3968       ENDDO
3969    ENDIF
3970
3971    IF ( nest_bound_n )  THEN
3972       j = ny + 1
3973       innor = -dx
3974       DO   i = nxl, nxr
3975          DO   k = nzb+1, nzt
3976             volume_flow_l(2) = volume_flow_l(2) + innor * v(k,j,i) * dzw(k)   &
3977                                     * MERGE( 1.0_wp, 0.0_wp,                  &
3978                                              BTEST( wall_flags_0(k,j,i), 2 ) )
3979          ENDDO
3980       ENDDO
3981    ENDIF
3982
3983#if defined( __parallel )
3984    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3985    CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL,         &
3986                        MPI_SUM, comm2d, ierr )
3987#else
3988    volume_flow(2) = volume_flow_l(2)
3989#endif
3990
3991!
3992!-- Sum up the volume flow through the top boundary
3993    volume_flow(3)   = 0.0_wp
3994    volume_flow_l(3) = 0.0_wp
3995    dxdy = dx * dy
3996    k = nzt
3997    DO   i = nxl, nxr
3998       DO   j = nys, nyn
3999          volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dxdy
4000       ENDDO
4001    ENDDO
4002
4003#if defined( __parallel )
4004    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
4005    CALL MPI_ALLREDUCE( volume_flow_l(3), volume_flow(3), 1, MPI_REAL,         &
4006                        MPI_SUM, comm2d, ierr )
4007#else
4008    volume_flow(3) = volume_flow_l(3)
4009#endif
4010
4011!
4012!-- Correct the top-boundary value of w
4013    w_lt = (volume_flow(1) + volume_flow(2) + volume_flow(3)) / area_t
4014    DO   i = nxl, nxr
4015       DO   j = nys, nyn
4016          DO  k = nzt, nzt + 1
4017             w(k,j,i) = w(k,j,i) + w_lt
4018          ENDDO
4019       ENDDO
4020    ENDDO
4021
4022 END SUBROUTINE pmci_ensure_nest_mass_conservation
4023
4024
4025
4026 SUBROUTINE pmci_synchronize
4027
4028#if defined( __parallel )
4029!
4030!-- Unify the time steps for each model and synchronize using
4031!-- MPI_ALLREDUCE with the MPI_MIN operator over all processes using
4032!-- the global communicator MPI_COMM_WORLD.
4033   
4034   IMPLICIT NONE
4035
4036   INTEGER(iwp)           :: ierr  !<
4037   REAL(wp), DIMENSION(1) :: dtl   !<
4038   REAL(wp), DIMENSION(1) :: dtg   !<
4039
4040   
4041   dtl(1) = dt_3d
4042   CALL MPI_ALLREDUCE( dtl, dtg, 1, MPI_REAL, MPI_MIN, MPI_COMM_WORLD, ierr )
4043   dt_3d  = dtg(1)
4044
4045#endif
4046 END SUBROUTINE pmci_synchronize
4047               
4048
4049
4050 SUBROUTINE pmci_set_swaplevel( swaplevel )
4051
4052!
4053!-- After each Runge-Kutta sub-timestep, alternately set buffer one or buffer
4054!-- two active
4055
4056    IMPLICIT NONE
4057
4058    INTEGER(iwp), INTENT(IN) ::  swaplevel  !< swaplevel (1 or 2) of PALM's
4059                                            !< timestep
4060
4061    INTEGER(iwp)            ::  child_id    !<
4062    INTEGER(iwp)            ::  m           !<
4063
4064#if defined( __parallel )
4065    DO  m = 1, SIZE( pmc_parent_for_child )-1
4066       child_id = pmc_parent_for_child(m)
4067       CALL pmc_s_set_active_data_array( child_id, swaplevel )
4068    ENDDO
4069#endif
4070 END SUBROUTINE pmci_set_swaplevel
4071
4072
4073
4074 SUBROUTINE pmci_datatrans( local_nesting_mode )
4075!
4076!-- This subroutine controls the nesting according to the nestpar
4077!-- parameter nesting_mode (two-way (default) or one-way) and the
4078!-- order of anterpolations according to the nestpar parameter
4079!-- nesting_datatransfer_mode (cascade, overlap or mixed (default)).
4080!-- Although nesting_mode is a variable of this model, pass it as
4081!-- an argument to allow for example to force one-way initialization
4082!-- phase.
4083
4084    IMPLICIT NONE
4085
4086    INTEGER(iwp)           ::  ierr   !<
4087    INTEGER(iwp)           ::  istat  !<
4088
4089    CHARACTER(LEN=*), INTENT(IN) ::  local_nesting_mode
4090
4091    IF ( TRIM( local_nesting_mode ) == 'one-way' )  THEN
4092
4093       CALL pmci_child_datatrans( parent_to_child )
4094       CALL pmci_parent_datatrans( parent_to_child )
4095
4096    ELSE
4097
4098       IF( nesting_datatransfer_mode == 'cascade' )  THEN
4099
4100          CALL pmci_child_datatrans( parent_to_child )
4101          CALL pmci_parent_datatrans( parent_to_child )
4102
4103          CALL pmci_parent_datatrans( child_to_parent )
4104          CALL pmci_child_datatrans( child_to_parent )
4105
4106       ELSEIF( nesting_datatransfer_mode == 'overlap')  THEN
4107
4108          CALL pmci_parent_datatrans( parent_to_child )
4109          CALL pmci_child_datatrans( parent_to_child )
4110
4111          CALL pmci_child_datatrans( child_to_parent )
4112          CALL pmci_parent_datatrans( child_to_parent )
4113
4114       ELSEIF( TRIM( nesting_datatransfer_mode ) == 'mixed' )  THEN
4115
4116          CALL pmci_parent_datatrans( parent_to_child )
4117          CALL pmci_child_datatrans( parent_to_child )
4118
4119          CALL pmci_parent_datatrans( child_to_parent )
4120          CALL pmci_child_datatrans( child_to_parent )
4121
4122       ENDIF
4123
4124    ENDIF
4125
4126 END SUBROUTINE pmci_datatrans
4127
4128
4129
4130 SUBROUTINE pmci_parent_datatrans( direction )
4131
4132    IMPLICIT NONE
4133
4134    INTEGER(iwp), INTENT(IN) ::  direction   !<
4135
4136#if defined( __parallel )
4137    INTEGER(iwp) ::  child_id    !<
4138    INTEGER(iwp) ::  i           !<
4139    INTEGER(iwp) ::  ierr        !<
4140    INTEGER(iwp) ::  j           !<
4141    INTEGER(iwp) ::  k           !<
4142    INTEGER(iwp) ::  m           !<
4143
4144    REAL(wp)               ::  waittime    !<
4145    REAL(wp), DIMENSION(1) ::  dtc         !<
4146    REAL(wp), DIMENSION(1) ::  dtl         !<
4147
4148
4149    DO  m = 1, SIZE( pmc_parent_for_child ) - 1
4150       child_id = pmc_parent_for_child(m)
4151       
4152       IF ( direction == parent_to_child )  THEN
4153          CALL cpu_log( log_point_s(71), 'pmc parent send', 'start' )
4154          CALL pmc_s_fillbuffer( child_id )
4155          CALL cpu_log( log_point_s(71), 'pmc parent send', 'stop' )
4156       ELSE
4157!
4158!--       Communication from child to parent
4159          CALL cpu_log( log_point_s(72), 'pmc parent recv', 'start' )
4160          child_id = pmc_parent_for_child(m)
4161          CALL pmc_s_getdata_from_buffer( child_id )
4162          CALL cpu_log( log_point_s(72), 'pmc parent recv', 'stop' )
4163!
4164!--       The anterpolated data is now available in u etc
4165          IF ( topography /= 'flat' )  THEN
4166!
4167!--          Inside buildings/topography reset velocities back to zero.
4168!--          Scalars (pt, q, s, km, kh, p, sa, ...) are ignored at
4169!--          present, maybe revise later.
4170!--          Resetting of e is removed as unnecessary since e is not
4171!--          anterpolated, and as incorrect since it overran the default
4172!--          Neumann condition (bc_e_b).
4173             DO   i = nxlg, nxrg
4174                DO   j = nysg, nyng
4175                   DO  k = nzb, nzt+1
4176                      u(k,j,i)  = MERGE( u(k,j,i), 0.0_wp,                     &
4177                                         BTEST( wall_flags_0(k,j,i), 1 ) )
4178                      v(k,j,i)  = MERGE( v(k,j,i), 0.0_wp,                     &
4179                                         BTEST( wall_flags_0(k,j,i), 2 ) )
4180                      w(k,j,i)  = MERGE( w(k,j,i), 0.0_wp,                     &
4181                                         BTEST( wall_flags_0(k,j,i), 3 ) )
4182!
4183!--                 TO_DO: zero setting of temperature within topography creates
4184!--                       wrong results
4185!                   pt(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
4186!                   IF ( humidity  .OR.  passive_scalar )  THEN
4187!                      q(nzb:nzb_s_inner(j,i),j,i) = 0.0_wp
4188!                   ENDIF
4189                   ENDDO
4190                ENDDO
4191             ENDDO
4192          ENDIF
4193       ENDIF
4194    ENDDO
4195
4196#endif
4197 END SUBROUTINE pmci_parent_datatrans
4198
4199
4200
4201 SUBROUTINE pmci_child_datatrans( direction )
4202
4203    IMPLICIT NONE
4204
4205    INTEGER(iwp), INTENT(IN) ::  direction   !<
4206
4207#if defined( __parallel )
4208    INTEGER(iwp) ::  ierr        !<
4209    INTEGER(iwp) ::  icl         !<
4210    INTEGER(iwp) ::  icr         !<
4211    INTEGER(iwp) ::  jcs         !<
4212    INTEGER(iwp) ::  jcn         !<
4213   
4214    REAL(wp), DIMENSION(1) ::  dtl         !<
4215    REAL(wp), DIMENSION(1) ::  dts         !<
4216
4217
4218    dtl = dt_3d
4219    IF ( cpl_id > 1 )  THEN
4220!
4221!--    Child domain boundaries in the parent indice space.
4222       icl = coarse_bound(1)
4223       icr = coarse_bound(2)
4224       jcs = coarse_bound(3)
4225       jcn = coarse_bound(4)
4226
4227       IF ( direction == parent_to_child )  THEN
4228
4229          CALL cpu_log( log_point_s(73), 'pmc child recv', 'start' )
4230          CALL pmc_c_getbuffer( )
4231          CALL cpu_log( log_point_s(73), 'pmc child recv', 'stop' )
4232
4233          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'start' )
4234          CALL pmci_interpolation
4235          CALL cpu_log( log_point_s(75), 'pmc interpolation', 'stop' )
4236
4237       ELSE
4238!
4239!--       direction == child_to_parent
4240          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'start' )
4241          CALL pmci_anterpolation
4242          CALL cpu_log( log_point_s(76), 'pmc anterpolation', 'stop' )
4243
4244          CALL cpu_log( log_point_s(74), 'pmc child send', 'start' )
4245          CALL pmc_c_putbuffer( )
4246          CALL cpu_log( log_point_s(74), 'pmc child send', 'stop' )
4247
4248       ENDIF
4249    ENDIF
4250
4251  CONTAINS
4252
4253   
4254    SUBROUTINE pmci_interpolation
4255
4256!
4257!--    A wrapper routine for all interpolation actions
4258     
4259       IMPLICIT NONE
4260
4261       INTEGER(iwp) ::  n          !< running index for number of chemical species
4262     
4263!
4264!--    In case of vertical nesting no interpolation is needed for the
4265!--    horizontal boundaries
4266       IF ( nesting_mode /= 'vertical' )  THEN
4267       
4268!
4269!--       Left border pe:
4270          IF ( nest_bound_l )  THEN
4271             
4272             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
4273                                       r1yo, r2yo, r1zo, r2zo,                 &
4274                                       logc_u_l, logc_ratio_u_l,               &
4275                                       logc_kbounds_u_l,                       &
4276                                       nzt_topo_nestbc_l, 'l', 'u' )
4277
4278             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
4279                                       r1yv, r2yv, r1zo, r2zo,                 &
4280                                       logc_v_l, logc_ratio_v_l,               &
4281                                       logc_kbounds_v_l,                       &               
4282                                       nzt_topo_nestbc_l, 'l', 'v' )
4283
4284             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
4285                                       r1yo, r2yo, r1zw, r2zw,                 &
4286                                       logc_w_l, logc_ratio_w_l,               &
4287                                       logc_kbounds_w_l,                       &
4288                                       nzt_topo_nestbc_l, 'l', 'w' )
4289
4290             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
4291                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
4292                     .NOT. constant_diffusion ) )  THEN
4293                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
4294                                          r1yo, r2yo, r1zo, r2zo,              &
4295                                          logc_w_l, logc_ratio_w_l,            &
4296                                          logc_kbounds_w_l,                    &
4297                                          nzt_topo_nestbc_l, 'l', 'e' )
4298             ENDIF
4299
4300             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
4301                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
4302                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
4303                                          logc_w_l, logc_ratio_w_l,            &
4304                                          logc_kbounds_w_l,                    &
4305                                          nzt_topo_nestbc_l, 'l', 's' )
4306             ENDIF
4307
4308             IF ( .NOT. neutral )  THEN
4309                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
4310                                          r1yo, r2yo, r1zo, r2zo,              &
4311                                          logc_w_l, logc_ratio_w_l,            &
4312                                          logc_kbounds_w_l,                    &               
4313                                          nzt_topo_nestbc_l, 'l', 's' )
4314             ENDIF
4315
4316             IF ( humidity )  THEN
4317
4318                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
4319                                          r1yo, r2yo, r1zo, r2zo,              &
4320                                          logc_w_l, logc_ratio_w_l,            &
4321                                          logc_kbounds_w_l,                    &
4322                                          nzt_topo_nestbc_l, 'l', 's' )
4323
4324                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4325                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
4326                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4327                                             logc_w_l, logc_ratio_w_l,         &
4328                                             logc_kbounds_w_l,                 &
4329                                             nzt_topo_nestbc_l, 'l', 's' ) 
4330
4331                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
4332                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4333                                             logc_w_l, logc_ratio_w_l,         &
4334                                             logc_kbounds_w_l,                 &
4335                                             nzt_topo_nestbc_l, 'l', 's' )         
4336                ENDIF
4337
4338                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4339                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
4340                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4341                                             logc_w_l, logc_ratio_w_l,         &
4342                                             logc_kbounds_w_l,                 &
4343                                             nzt_topo_nestbc_l, 'l', 's' ) 
4344
4345                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
4346                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4347                                             logc_w_l, logc_ratio_w_l,         &
4348                                             logc_kbounds_w_l,                 &               
4349                                             nzt_topo_nestbc_l, 'l', 's' )             
4350                ENDIF
4351
4352             ENDIF
4353
4354             IF ( passive_scalar )  THEN
4355                CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
4356                                          r1yo, r2yo, r1zo, r2zo,              &
4357                                          logc_w_l, logc_ratio_w_l,            &
4358                                          logc_kbounds_w_l,                    & 
4359                                          nzt_topo_nestbc_l, 'l', 's' )
4360             ENDIF
4361
4362             IF ( air_chemistry )  THEN
4363                DO  n = 1, nspec
4364                   CALL pmci_interp_tril_lr( chem_species(n)%conc,             &
4365                                             chem_spec_c(:,:,:,n),             &
4366                                             ico, jco, kco, r1xo, r2xo,        &
4367                                             r1yo, r2yo, r1zo, r2zo,           &
4368                                             logc_w_l, logc_ratio_w_l,         &
4369                                             logc_kbounds_w_l,                 & 
4370                                             nzt_topo_nestbc_l, 'l', 's' )
4371                ENDDO
4372             ENDIF
4373
4374          ENDIF
4375!
4376!--       Right border pe
4377          IF ( nest_bound_r )  THEN
4378             
4379             CALL pmci_interp_tril_lr( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
4380                                       r1yo, r2yo, r1zo, r2zo,                 &
4381                                       logc_u_r, logc_ratio_u_r,               &
4382                                       logc_kbounds_u_r,                       &
4383                                       nzt_topo_nestbc_r, 'r', 'u' )
4384
4385             CALL pmci_interp_tril_lr( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
4386                                       r1yv, r2yv, r1zo, r2zo,                 &
4387                                       logc_v_r, logc_ratio_v_r,               &
4388                                       logc_kbounds_v_r,                       &
4389                                       nzt_topo_nestbc_r, 'r', 'v' )
4390
4391             CALL pmci_interp_tril_lr( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
4392                                       r1yo, r2yo, r1zw, r2zw,                 &
4393                                       logc_w_r, logc_ratio_w_r,               &
4394                                       logc_kbounds_w_r,                       &
4395                                       nzt_topo_nestbc_r, 'r', 'w' )
4396
4397             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
4398                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
4399                     .NOT. constant_diffusion ) )  THEN
4400                CALL pmci_interp_tril_lr( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
4401                                          r1yo,r2yo, r1zo, r2zo,               &
4402                                          logc_w_r, logc_ratio_w_r,            &
4403                                          logc_kbounds_w_r,                    &
4404                                          nzt_topo_nestbc_r, 'r', 'e' )
4405
4406             ENDIF
4407
4408             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
4409                CALL pmci_interp_tril_lr( diss,  dissc,  ico, jco, kco, r1xo,  &
4410                                          r2xo, r1yo,r2yo, r1zo, r2zo,         &
4411                                          logc_w_r, logc_ratio_w_r,            &
4412                                          logc_kbounds_w_r,                    &
4413                                          nzt_topo_nestbc_r, 'r', 's' )
4414
4415             ENDIF
4416
4417             IF ( .NOT. neutral )  THEN
4418                CALL pmci_interp_tril_lr( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
4419                                          r1yo, r2yo, r1zo, r2zo,              &
4420                                          logc_w_r, logc_ratio_w_r,            &
4421                                          logc_kbounds_w_r,                    &
4422                                          nzt_topo_nestbc_r, 'r', 's' )
4423             ENDIF
4424
4425             IF ( humidity )  THEN
4426                CALL pmci_interp_tril_lr( q, q_c, ico, jco, kco, r1xo, r2xo,   &
4427                                          r1yo, r2yo, r1zo, r2zo,              &
4428                                          logc_w_r, logc_ratio_w_r,            &
4429                                          logc_kbounds_w_r,                    &
4430                                          nzt_topo_nestbc_r, 'r', 's' )
4431
4432                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4433
4434                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
4435                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4436                                             logc_w_r, logc_ratio_w_r,         &
4437                                             logc_kbounds_w_r,                 &
4438                                             nzt_topo_nestbc_r, 'r', 's' ) 
4439     
4440                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
4441                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4442                                             logc_w_r, logc_ratio_w_r,         &
4443                                             logc_kbounds_w_r,                 &
4444                                             nzt_topo_nestbc_r, 'r', 's' )
4445
4446
4447                ENDIF
4448
4449                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4450
4451     
4452                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
4453                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4454                                             logc_w_r, logc_ratio_w_r,         &
4455                                             logc_kbounds_w_r,                 &
4456                                             nzt_topo_nestbc_r, 'r', 's' ) 
4457
4458                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
4459                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4460                                             logc_w_r, logc_ratio_w_r,         &
4461                                             logc_kbounds_w_r,                 &
4462                                             nzt_topo_nestbc_r, 'r', 's' ) 
4463
4464                ENDIF
4465
4466             ENDIF
4467
4468             IF ( passive_scalar )  THEN
4469                CALL pmci_interp_tril_lr( s, sc, ico, jco, kco, r1xo, r2xo,    &
4470                                          r1yo, r2yo, r1zo, r2zo,              &
4471                                          logc_w_r, logc_ratio_w_r,            &
4472                                          logc_kbounds_w_r,                    &
4473                                          nzt_topo_nestbc_r, 'r', 's' )
4474             ENDIF
4475
4476             IF ( air_chemistry )  THEN
4477                DO  n = 1, nspec
4478                   CALL pmci_interp_tril_lr( chem_species(n)%conc,             &
4479                                             chem_spec_c(:,:,:,n),             &
4480                                             ico, jco, kco, r1xo, r2xo,        &
4481                                             r1yo, r2yo, r1zo, r2zo,           &
4482                                             logc_w_r, logc_ratio_w_r,         &
4483                                             logc_kbounds_w_r,                 &
4484                                             nzt_topo_nestbc_r, 'r', 's' )
4485                ENDDO
4486             ENDIF
4487          ENDIF
4488!
4489!--       South border pe
4490          IF ( nest_bound_s )  THEN
4491
4492             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
4493                                       r1yo, r2yo, r1zo, r2zo,                 &
4494                                       logc_u_s, logc_ratio_u_s,               &
4495                                       logc_kbounds_u_s,                       &
4496                                       nzt_topo_nestbc_s, 's', 'u' )
4497
4498             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
4499                                       r1yv, r2yv, r1zo, r2zo,                 &
4500                                       logc_v_s, logc_ratio_v_s,               &
4501                                       logc_kbounds_v_s,                       &
4502                                       nzt_topo_nestbc_s, 's', 'v' )
4503
4504             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
4505                                       r1yo, r2yo, r1zw, r2zw,                 &
4506                                       logc_w_s, logc_ratio_w_s,               &
4507                                       logc_kbounds_w_s,                       &
4508                                       nzt_topo_nestbc_s, 's','w' )
4509
4510             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   &
4511                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
4512                     .NOT. constant_diffusion ) )  THEN
4513                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
4514                                          r1yo, r2yo, r1zo, r2zo,              &
4515                                          logc_w_s, logc_ratio_w_s,            &
4516                                          logc_kbounds_w_s,                    &
4517                                          nzt_topo_nestbc_s, 's', 'e' )
4518
4519             ENDIF
4520
4521             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
4522                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
4523                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
4524                                          logc_w_s, logc_ratio_w_s,            &
4525                                          logc_kbounds_w_s,                    &
4526                                          nzt_topo_nestbc_s, 's', 's' )
4527
4528             ENDIF
4529
4530             IF ( .NOT. neutral )  THEN
4531                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
4532                                          r1yo, r2yo, r1zo, r2zo,              &
4533                                          logc_w_s, logc_ratio_w_s,            &
4534                                          logc_kbounds_w_s,                    &
4535                                          nzt_topo_nestbc_s, 's', 's' )
4536             ENDIF
4537
4538             IF ( humidity )  THEN
4539                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
4540                                          r1yo,r2yo, r1zo, r2zo,               &
4541                                          logc_w_s, logc_ratio_w_s,            &
4542                                          logc_kbounds_w_s,                    &
4543                                          nzt_topo_nestbc_s, 's', 's' )
4544
4545                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4546
4547                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
4548                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
4549                                             logc_w_s, logc_ratio_w_s,         &
4550                                             logc_kbounds_w_s,                 &
4551                                             nzt_topo_nestbc_s, 's', 's' )
4552
4553                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
4554                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
4555                                             logc_w_s, logc_ratio_w_s,         &
4556                                             logc_kbounds_w_s,                 &
4557                                             nzt_topo_nestbc_s, 's', 's' )
4558
4559                ENDIF
4560
4561                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4562
4563                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
4564                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
4565                                             logc_w_s, logc_ratio_w_s,         &
4566                                             logc_kbounds_w_s,                 &
4567                                             nzt_topo_nestbc_s, 's', 's' )
4568
4569                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
4570                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
4571                                             logc_w_s, logc_ratio_w_s,         &
4572                                             logc_kbounds_w_s,                 &
4573                                             nzt_topo_nestbc_s, 's', 's' )
4574
4575                ENDIF
4576
4577             ENDIF
4578
4579             IF ( passive_scalar )  THEN
4580                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
4581                                          r1yo,r2yo, r1zo, r2zo,               &
4582                                          logc_w_s, logc_ratio_w_s,            &
4583                                          logc_kbounds_w_s,                    &
4584                                          nzt_topo_nestbc_s, 's', 's' )
4585             ENDIF
4586
4587             IF ( air_chemistry )  THEN
4588                DO  n = 1, nspec
4589                   CALL pmci_interp_tril_sn( chem_species(n)%conc,             &
4590                                             chem_spec_c(:,:,:,n),             &
4591                                             ico, jco, kco, r1xo, r2xo,        &
4592                                             r1yo, r2yo, r1zo, r2zo,           &
4593                                             logc_w_s, logc_ratio_w_s,         &
4594                                             logc_kbounds_w_s,                 &
4595                                             nzt_topo_nestbc_s, 's', 's' )
4596                ENDDO
4597             ENDIF
4598          ENDIF
4599!
4600!--       North border pe
4601          IF ( nest_bound_n )  THEN
4602             
4603             CALL pmci_interp_tril_sn( u,  uc,  icu, jco, kco, r1xu, r2xu,     &
4604                                       r1yo, r2yo, r1zo, r2zo,                 &
4605                                       logc_u_n, logc_ratio_u_n,               &
4606                                       logc_kbounds_u_n,                       &
4607                                       nzt_topo_nestbc_n, 'n', 'u' )
4608
4609             CALL pmci_interp_tril_sn( v,  vc,  ico, jcv, kco, r1xo, r2xo,     &
4610                                       r1yv, r2yv, r1zo, r2zo,                 &
4611                                       logc_v_n, logc_ratio_v_n,               &
4612                                       logc_kbounds_v_n,                       &
4613                                       nzt_topo_nestbc_n, 'n', 'v' )
4614
4615             CALL pmci_interp_tril_sn( w,  wc,  ico, jco, kcw, r1xo, r2xo,     &
4616                                       r1yo, r2yo, r1zw, r2zw,                 &
4617                                       logc_w_n, logc_ratio_w_n,               &
4618                                       logc_kbounds_w_n,                       &
4619                                       nzt_topo_nestbc_n, 'n', 'w' )
4620
4621             IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.   & 
4622                  (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.    &
4623                     .NOT. constant_diffusion ) )  THEN
4624                CALL pmci_interp_tril_sn( e,  ec,  ico, jco, kco, r1xo, r2xo,  &
4625                                          r1yo, r2yo, r1zo, r2zo,              &
4626                                          logc_w_n, logc_ratio_w_n,            &
4627                                          logc_kbounds_w_n,                    &
4628                                          nzt_topo_nestbc_n, 'n', 'e' )
4629
4630             ENDIF
4631
4632             IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
4633                CALL pmci_interp_tril_sn( diss, dissc,  ico, jco, kco, r1xo,   &
4634                                          r2xo, r1yo, r2yo, r1zo, r2zo,        &
4635                                          logc_w_n, logc_ratio_w_n,            &
4636                                          logc_kbounds_w_n,                    &
4637                                          nzt_topo_nestbc_n, 'n', 's' )
4638
4639             ENDIF
4640
4641             IF ( .NOT. neutral )  THEN
4642                CALL pmci_interp_tril_sn( pt, ptc, ico, jco, kco, r1xo, r2xo,  &
4643                                          r1yo, r2yo, r1zo, r2zo,              &
4644                                          logc_w_n, logc_ratio_w_n,            &
4645                                          logc_kbounds_w_n,                    &
4646                                          nzt_topo_nestbc_n, 'n', 's' )
4647             ENDIF
4648
4649             IF ( humidity )  THEN
4650                CALL pmci_interp_tril_sn( q, q_c, ico, jco, kco, r1xo, r2xo,   &
4651                                          r1yo, r2yo, r1zo, r2zo,              &
4652                                          logc_w_n, logc_ratio_w_n,            &
4653                                          logc_kbounds_w_n,                    &
4654                                          nzt_topo_nestbc_n, 'n', 's' )
4655
4656                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4657
4658                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
4659                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4660                                             logc_w_n, logc_ratio_w_n,         &
4661                                             logc_kbounds_w_n,                 &
4662                                             nzt_topo_nestbc_n, 'n', 's' )
4663
4664                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
4665                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4666                                             logc_u_n, logc_ratio_u_n,         &
4667                                             logc_kbounds_w_n,                 &
4668                                             nzt_topo_nestbc_n, 'n', 's' )
4669
4670                ENDIF
4671
4672                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4673
4674                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
4675                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4676                                             logc_w_n, logc_ratio_w_n,         &
4677                                             logc_kbounds_w_n,                 &
4678                                             nzt_topo_nestbc_n, 'n', 's' )
4679
4680                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
4681                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
4682                                             logc_w_n, logc_ratio_w_n,         &
4683                                             logc_kbounds_w_n,                 &
4684                                             nzt_topo_nestbc_n, 'n', 's' )
4685
4686                ENDIF
4687
4688             ENDIF
4689
4690             IF ( passive_scalar )  THEN
4691                CALL pmci_interp_tril_sn( s, sc, ico, jco, kco, r1xo, r2xo,    &
4692                                          r1yo, r2yo, r1zo, r2zo,              &
4693                                          logc_w_n, logc_ratio_w_n,            &
4694                                          logc_kbounds_w_n,                    &
4695                                          nzt_topo_nestbc_n, 'n', 's' )
4696             ENDIF
4697
4698             IF ( air_chemistry )  THEN
4699                DO  n = 1, nspec
4700                   CALL pmci_interp_tril_sn( chem_species(n)%conc,             &
4701                                             chem_spec_c(:,:,:,n),             &
4702                                             ico, jco, kco, r1xo, r2xo,        &
4703                                             r1yo, r2yo, r1zo, r2zo,           &
4704                                             logc_w_n, logc_ratio_w_n,         &
4705                                             logc_kbounds_w_n,                 &
4706                                             nzt_topo_nestbc_n, 'n', 's' )
4707                ENDDO
4708             ENDIF
4709          ENDIF
4710
4711       ENDIF       ! IF ( nesting_mode /= 'vertical' )
4712!
4713!--    All PEs are top-border PEs
4714       CALL pmci_interp_tril_t( u,  uc,  icu, jco, kco, r1xu, r2xu, r1yo,      &
4715                                r2yo, r1zo, r2zo, 'u' )
4716       CALL pmci_interp_tril_t( v,  vc,  ico, jcv, kco, r1xo, r2xo, r1yv,      &
4717                                r2yv, r1zo, r2zo, 'v' )
4718       CALL pmci_interp_tril_t( w,  wc,  ico, jco, kcw, r1xo, r2xo, r1yo,      &
4719                                r2yo, r1zw, r2zw, 'w' )
4720
4721       IF ( (        rans_mode_parent  .AND.         rans_mode )  .OR.         &
4722            (  .NOT. rans_mode_parent  .AND.  .NOT.  rans_mode  .AND.          &
4723               .NOT. constant_diffusion ) )  THEN
4724          CALL pmci_interp_tril_t( e,  ec,  ico, jco, kco, r1xo, r2xo, r1yo,   &
4725                                   r2yo, r1zo, r2zo, 'e' )
4726       ENDIF
4727
4728       IF ( rans_mode_parent  .AND.  rans_mode  .AND.  rans_tke_e )  THEN
4729          CALL pmci_interp_tril_t( diss, dissc, ico, jco, kco, r1xo, r2xo,     &
4730                                   r1yo, r2yo, r1zo, r2zo, 's' )
4731       ENDIF
4732
4733       IF ( .NOT. neutral )  THEN
4734          CALL pmci_interp_tril_t( pt, ptc, ico, jco, kco, r1xo, r2xo, r1yo,   &
4735                                   r2yo, r1zo, r2zo, 's' )
4736       ENDIF
4737
4738       IF ( humidity )  THEN
4739
4740          CALL pmci_interp_tril_t( q, q_c, ico, jco, kco, r1xo, r2xo, r1yo,    &
4741                                   r2yo, r1zo, r2zo, 's' )
4742
4743          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4744
4745             CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
4746                                      r2yo, r1zo, r2zo, 's' )
4747
4748             CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
4749                                      r2yo, r1zo, r2zo, 's' )
4750
4751          ENDIF
4752
4753          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4754
4755
4756             CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
4757                                      r2yo, r1zo, r2zo, 's' )
4758
4759             CALL pmci_interp_tril_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,&
4760                                      r2yo, r1zo, r2zo, 's' )
4761
4762          ENDIF
4763
4764       ENDIF
4765
4766       IF ( passive_scalar )  THEN
4767          CALL pmci_interp_tril_t( s, sc, ico, jco, kco, r1xo, r2xo, r1yo,     &
4768                                   r2yo, r1zo, r2zo, 's' )
4769       ENDIF
4770
4771       IF ( air_chemistry )  THEN
4772          DO  n = 1, nspec
4773             CALL pmci_interp_tril_t( chem_species(n)%conc,                    &
4774                                      chem_spec_c(:,:,:,n),                    &
4775                                      ico, jco, kco, r1xo, r2xo,               &
4776                                      r1yo, r2yo, r1zo, r2zo,                  &
4777                                      's' )
4778          ENDDO
4779       ENDIF
4780
4781   END SUBROUTINE pmci_interpolation
4782
4783
4784
4785   SUBROUTINE pmci_anterpolation
4786
4787!
4788!--   A wrapper routine for all anterpolation actions.
4789!--   Note that TKE is not anterpolated.
4790      IMPLICIT NONE
4791
4792      INTEGER(iwp) ::  n          !< running index for number of chemical species
4793
4794
4795
4796      CALL pmci_anterp_tophat( u,  uc,  kctu, iflu, ifuu, jflo, jfuo, kflo,    &
4797                               kfuo, ijfc_u, kfc_s, 'u' )
4798      CALL pmci_anterp_tophat( v,  vc,  kctu, iflo, ifuo, jflv, jfuv, kflo,    &
4799                               kfuo, ijfc_v, kfc_s, 'v' )
4800      CALL pmci_anterp_tophat( w,  wc,  kctw, iflo, ifuo, jflo, jfuo, kflw,    &
4801                               kfuw, ijfc_s, kfc_w, 'w' )
4802!
4803!--   Anterpolation of TKE and dissipation rate if parent and child are in
4804!--   RANS mode.
4805      IF ( rans_mode_parent  .AND.  rans_mode )  THEN
4806         CALL pmci_anterp_tophat( e, ec, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
4807                                  kfuo, ijfc_s, kfc_s, 'e' )
4808!
4809!--      Anterpolation of dissipation rate only if TKE-e closure is applied.
4810         IF ( rans_tke_e )  THEN
4811            CALL pmci_anterp_tophat( diss, dissc, kctu, iflo, ifuo, jflo, jfuo,&
4812                                     kflo, kfuo, ijfc_s, kfc_s, 'diss' )
4813         ENDIF
4814
4815      ENDIF
4816
4817      IF ( .NOT. neutral )  THEN
4818         CALL pmci_anterp_tophat( pt, ptc, kctu, iflo, ifuo, jflo, jfuo, kflo, &
4819                                  kfuo, ijfc_s, kfc_s, 'pt' )
4820      ENDIF
4821
4822      IF ( humidity )  THEN
4823
4824         CALL pmci_anterp_tophat( q, q_c, kctu, iflo, ifuo, jflo, jfuo, kflo,  &
4825                                  kfuo, ijfc_s, kfc_s, 'q' )
4826
4827         IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
4828
4829            CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
4830                                     kflo, kfuo, ijfc_s, kfc_s, 'qc' )
4831
4832            CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
4833                                     kflo, kfuo, ijfc_s, kfc_s,  'nc' )
4834
4835         ENDIF
4836
4837         IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
4838
4839            CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
4840                                     kflo, kfuo, ijfc_s, kfc_s, 'qr' )
4841
4842            CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &
4843                                     kflo, kfuo, ijfc_s, kfc_s, 'nr' )
4844
4845         ENDIF
4846
4847      ENDIF
4848
4849      IF ( passive_scalar )  THEN
4850         CALL pmci_anterp_tophat( s, sc, kctu, iflo, ifuo, jflo, jfuo, kflo,   &
4851                                  kfuo, ijfc_s, kfc_s, 's' )
4852      ENDIF
4853
4854      IF ( air_chemistry )  THEN
4855         DO  n = 1, nspec
4856            CALL pmci_anterp_tophat( chem_species(n)%conc,                     &
4857                                     chem_spec_c(:,:,:,n),                     &
4858                                     kctu, iflo, ifuo, jflo, jfuo, kflo,       &
4859                                     kfuo, ijfc_s, kfc_s, 's' )
4860         ENDDO
4861      ENDIF
4862
4863   END SUBROUTINE pmci_anterpolation
4864
4865
4866
4867   SUBROUTINE pmci_interp_tril_lr( f, fc, ic, jc, kc, r1x, r2x, r1y, r2y, r1z, &
4868                                   r2z, logc, logc_ratio, logc_kbounds,        &
4869                                   nzt_topo_nestbc, edge, var )
4870!
4871!--   Interpolation of ghost-node values used as the child-domain boundary
4872!--   conditions. This subroutine handles the left and right boundaries. It is
4873!--   based on trilinear interpolation.
4874
4875      IMPLICIT NONE
4876
4877      REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                      &
4878                                      INTENT(INOUT) ::  f       !<
4879      REAL(wp), DIMENSION(0:cg%nz+1,jcs:jcn,icl:icr),                          &
4880                                      INTENT(IN)    ::  fc      !<
4881      REAL(wp), DIMENSION(1:2,0:ncorr-1,nzb:nzt_topo_nestbc,nys:nyn),          &
4882                                      INTENT(IN)    ::  logc_ratio   !<
4883      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r1x     !<
4884      REAL(wp), DIMENSION(nxlg:nxrg), INTENT(IN)    ::  r2x     !<
4885      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r1y     !<
4886      REAL(wp), DIMENSION(nysg:nyng), INTENT(IN)    ::  r2y     !<
4887      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r1z     !<
4888      REAL(wp), DIMENSION(nzb:nzt+1), INTENT(IN)    ::  r2z     !<
4889     
4890      INTEGER(iwp), DIMENSION(nxlg:nxrg), INTENT(IN)           ::  ic     !<
4891      INTEGER(iwp), DIMENSION(nysg:nyng), INTENT(IN)           ::  jc     !<
4892      INTEGER(iwp), DIMENSION(nzb:nzt+1), INTENT(IN)           ::  kc     !<
4893      INTEGER(iwp), DIMENSION(1:2,nzb:nzt_topo_nestbc,nys:nyn),                &
4894                                          INTENT(IN)           ::  logc   !<
4895      INTEGER(iwp), DIMENSION(1:2,nys:nyn), INTENT(IN)         ::  logc_kbounds !<
4896      INTEGER(iwp) ::  nzt_topo_nestbc   !<
4897
4898      CHARACTER(LEN=1), INTENT(IN) ::  edge   !<
4899      CHARACTER(LEN=1), INTENT(IN) ::  var    !<
4900
4901      INTEGER(iwp) ::  i        !<
4902      INTEGER(iwp) ::  ib       !<
4903      INTEGER(iwp) ::  ibgp     !<
4904      INTEGER(iwp) ::  iw       !<
4905      INTEGER(iwp) ::  j        !<
4906      INTEGER(iwp) ::  jco      !<
4907      INTEGER(iwp) ::  jcorr    !<
4908      INTEGER(iwp) ::  jinc     !<
4909      INTEGER(iwp) ::  jw       !<
4910      INTEGER(iwp) ::  j1       !<
4911      INTEGER(iwp) ::  k        !<
4912      INTEGER(iwp) ::  k_wall   !< vertical index of topography top
4913      INTEGER(iwp) ::  kco      !<
4914      INTEGER(iwp) ::  kcorr    !<
4915      INTEGER(iwp) ::  k1       !<
4916      INTEGER(iwp) ::  l        !<
4917      INTEGER(iwp) ::  m        !<
4918      INTEGER(iwp) ::  n        !<
4919      INTEGER(iwp) ::  kbc      !<
4920     
4921      REAL(wp) ::  coarse_dx   !<
4922      REAL(wp) ::  coarse_dy   !<
4923      REAL(wp) ::  coarse_dz   !<
4924      REAL(wp) ::  fkj         !<
4925      REAL(wp) ::  fkjp        !<
4926      REAL(wp) ::  fkpj        !<
4927      REAL(wp) ::  fkpjp       !<
4928      REAL(wp) ::  fk          !<
4929      REAL(wp) ::  fkp         !<
4930     
4931!
4932!--   Check which edge is to be handled
4933      IF ( edge == 'l' )  THEN
4934!
4935!--      For u, nxl is a ghost node, but not for the other variables
4936         IF ( var == 'u' )  THEN
4937            = nxl
4938            ib = nxl - 1 
4939         ELSE
4940            = nxl - 1
4941            ib = nxl - 2
4942         ENDIF
4943      ELSEIF ( edge == 'r' )  THEN
4944         = nxr + 1
4945         ib = nxr + 2
4946      ENDIF
4947     
4948      DO  j = nys, nyn+1
4949!
4950!--      Determine vertical index of topography top at grid point (j,i)
4951         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
4952
4953         DO  k = k_wall, nzt+1
4954            l = ic(i)
4955            m = jc(j)
4956            n = kc(k)
4957            fkj      = r1x(i) * fc(n,m,l)     + r2x(i) * fc(n,m,l+1)
4958            fkjp     = r1x(i) * fc(n,m+1,l)   + r2x(i) * fc(n,m+1,l+1)
4959            fkpj     = r1x(i) * fc(n+1,m,l)   + r2x(i) * fc(n+1,m,l+1)
4960            fkpjp    = r1x(i) * fc(n+1,m+1,l) + r2x(i) * fc(n+1,m+1,l+1)
4961            fk       = r1y(j) * fkj  + r2y(j) * fkjp
4962            fkp      = r1y(j) * fkpj + r2y(j) * fkpjp
4963            f(k,j,i) = r1z(k) * fk   + r2z(k) * fkp
4964         ENDDO
4965      ENDDO
4966!
4967!--   Generalized log-law-correction algorithm.
4968!--   Doubly two-dimensional index arrays logc(1:2,:,:) and log-ratio arrays
4969!--   logc_ratio(1:2,0:ncorr-1,:,:) have been precomputed in subroutine
4970!--   pmci_init_loglaw_correction.
4971!
4972!--   Solid surface below the node
4973      IF ( constant_flux_layer .AND. ( var == 'u' .OR. var == 'v' ) )  THEN
4974         DO  j = nys, nyn
4975!
4976!--         Determine vertical index of topography top at grid point (j,i)
4977            k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
4978
4979            k = k_wall+1
4980            IF ( ( logc(1,k,j) /= 0 )  .AND.  ( logc(2,k,j) == 0 ) )  THEN
4981               k1 = logc(1,k,j)
4982               DO  kcorr = 0, ncorr - 1
4983                  kco = k + kcorr
4984                  f(kco,j,i) = logc_ratio(1,kcorr,k,j) * f(k1,j,i)
4985               ENDDO
4986            ENDIF
4987         ENDDO
4988      ENDIF
4989!
4990!--   In case of non-flat topography, also vertical walls and corners need to be
4991!--   treated. Only single and double wall nodes are corrected. Triple and
4992!--   higher-multiple wall nodes are not corrected as the log law would not be
4993!--   valid anyway in such locations.
4994      IF ( topography /= 'flat' )  THEN
4995
4996         IF ( constant_flux_layer .AND. ( var == 'u' .OR. var == 'w' ) )  THEN           
4997!
4998!--         Solid surface only on south/north side of the node                   
4999            DO  j = nys, nyn
5000               DO  k = logc_kbounds(1,j), logc_kbounds(2,j)   
5001                  IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) == 0 ) )  THEN
5002!
5003!--                  Direction of the wall-normal index is carried in as the
5004!--                  sign of logc
5005                     jinc = SIGN( 1, logc(2,k,j) )
5006                     j1   = ABS( logc(2,k,j) )
5007                     DO  jcorr = 0, ncorr-1
5008                        jco = j + jinc * jcorr
5009                        IF ( jco >= nys .AND. jco <= nyn )  THEN
5010                           f(k,jco,i) = logc_ratio(2,jcorr,k,j) * f(k,j1,i)
5011                        ENDIF
5012                     ENDDO
5013                  ENDIF
5014               ENDDO
5015            ENDDO
5016         ENDIF
5017!
5018!--      Solid surface on both below and on south/north side of the node           
5019         IF ( constant_flux_layer .AND. var == 'u' )  THEN
5020            DO  j = nys, nyn
5021               k = logc_kbounds(1,j)
5022               IF ( ( logc(2,k,j) /= 0 )  .AND.  ( logc(1,k,j) /= 0 ) )  THEN
5023                  k1   = logc(1,k,j)                 
5024                  jinc = SIGN( 1, logc(2,k,j) )
5025                  j1   = ABS( logc(2,k,j) )                 
5026                  DO  jcorr = 0, ncorr-1
5027                     jco = j + jinc * jcorr
5028                     IF ( jco >= nys .AND. jco <= nyn )  THEN
5029                        DO  kcorr = 0, ncorr-1
5030                           kco = k + kcorr
5031                           f(kco,jco,i) = 0.5_wp * ( logc_ratio(1,kcorr,k,j) * &
5032                                                     f(k1,j,i)                 &
5033                                                   + logc_ratio(2,jcorr,k,j) * &
5034                                                     f(k,j1,i) )
5035                        ENDDO
5036                     ENDIF
5037                  ENDDO
5038               ENDIF
5039            ENDDO
5040         ENDIF
5041
5042      ENDIF  ! ( topography /= 'flat' )
5043!
5044!--   Rescale if f is the TKE.
5045      IF ( var == 'e')  THEN
5046         IF ( edge == 'l' )  THEN
5047            DO  j = nys, nyn + 1
5048!
5049!--            Determine vertical index of topography top at grid point (j,i)
5050               k_wall = get_topography_top_index_ji( j, i, 's' )
5051
5052               DO  k = k_wall, nzt + 1
5053                  f(k,j,i) = tkefactor_l(k,j) * f(k,j,i)
5054               ENDDO
5055            ENDDO
5056         ELSEIF ( edge == 'r' )  THEN           
5057            DO  j = nys, nyn+1
5058!
5059!--            Determine vertical index of topography top at grid point (j,i)
5060               k_wall = get_topography_top_index_ji( j, i, 's' )
5061
5062               DO  k = k_wall, nzt+1
5063                  f(k,j,i) = tkefactor_r(k,j) * f(k,j,i)
5064               ENDDO
5065            ENDDO
5066         ENDIF
5067      ENDIF
5068!
5069!--   Store the boundary values also into the other redundant ghost node layers.
5070!--   Please note, in case of only one ghost node layer, e.g. for the PW
5071!--   scheme, the following loops will not be entered.
5072      IF ( edge == 'l' )  THEN
5073         DO  ibgp = -nbgp, ib
5074            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
5075         ENDDO
5076      ELSEIF ( edge == 'r' )  THEN
5077         DO  ibgp = ib, nx+nbgp
5078            f(0:nzt+1,nysg:nyng,ibgp) = f(0:nzt+1,nysg:nyng,i)
5079         ENDDO
5080      ENDIF
5081
5082<