Home

Context Navigation

pmc_interface_mod.f90 @ 2938

Last change on this file since 2938 was 2938, checked in by suehring, 7 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/branches/palm4u/SOURCE/pmc_interface_mod.f90	2540-2692
/palm/trunk/SOURCE/pmc_interface_mod.f90	merged	eligible
/palm/branches/fricke/SOURCE/pmc_interface_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90	296-409
/palm/branches/suehring/pmc_interface_mod.f90	423-666

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

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

Context Navigation

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

Download in other formats: