Home

Context Navigation

pmc_interface_mod.f90 @ 2795

Last change on this file since 2795 was 2795, checked in by hellstea, 5 years ago

Bugfix in nesting anterpolation relaxation functions

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.6 KB

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