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

Last change on this file since 3021 was 3021, checked in by maronga, 3 years ago

bugfixes for nested runs

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