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

Last change on this file since 2938 was 2938, checked in by suehring, 4 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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