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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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