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

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

Revision history corrected

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