Home

Context Navigation

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

Last change on this file since 3979 was 3979, checked in by hellstea, 6 years ago
Bugfix in nesting interpolation pmci_interp_1sto_sn
Property svn:keywords set to `Id`
File size: 210.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |