Home

Context Navigation

pmc_interface_mod.f90 @ 2773

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

Nesting for chemical species implemented; Bugfix passive scalar boundary condition after anterpolation; Timeseries output of surface temperature; Enable initialization of 3D topography (was commented out so far)

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