Home

Context Navigation

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

Last change on this file since 3883 was 3883, checked in by hellstea, 5 years ago
Checks and error messages improved and extended. Number of variables renamed etc
Property svn:keywords set to `Id`
File size: 200.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |