Home

Context Navigation

pmc_interface_mod.f90 @ 3984

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