Home

Context Navigation

pmc_interface_mod.f90 @ 3592

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