Home

Context Navigation

pmc_interface_mod.f90 @ 3215

Last change on this file since 3215 was 3215, checked in by suehring, 5 years ago
changes for commit 3209 documented
Property svn:keywords set to `Id`
File size: 234.2 KB

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