Home

Context Navigation

pmc_interface_mod.f90 @ 3885

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