Home

Context Navigation

pmc_interface_mod.f90 @ 3182

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