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

Last change on this file since 3022 was 3022, checked in by suehring, 3 years ago

Revise recent bugfix in nested runs at left and south boundary; bugfix in advection of u in case of OpenMP parallelization; bugfix in plant transpiration

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