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

Last change on this file since 3655 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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