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

Last change on this file since 3768 was 3741, checked in by hellstea, 2 years ago

Bugfix in nest initialization and interpolations

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