Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |