Home

Context Navigation

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

Last change on this file since 4682 was 4665, checked in by hellstea, 5 years ago
Interpolation and anterpolation subroutines renamed and all missing subroutine description comments added (pmc_interface_mod). Format descriptor 102 slightly modified (cpulog_mod).
Property svn:keywords set to `Id`
File size: 234.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |