Home

Context Navigation

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

Last change on this file since 3876 was 3876, checked in by knoop, 5 years ago
Moved "photolysis_scheme", "chem_species" and "phot_frequen" to chem_modules
Property svn:keywords set to `Id`
File size: 187.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |