Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |