Home

Context Navigation

pmc_interface_mod.f90 @ 3976

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