Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |