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

Last change on this file since 3241 was 3241, checked in by raasch, 3 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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