Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |