Home

Context Navigation

pmc_interface_mod.f90 @ 3697

Last change on this file since 3697 was 3697, checked in by hellstea, 4 years ago
bugfix in child initialization in pmc_interface_mod
Property svn:keywords set to `Id`
File size: 343.8 KB

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