Home

Context Navigation

pmc_interface_mod.f90 @ 4260

Last change on this file since 4260 was 4260, checked in by hellstea, 5 years ago
Rest of the grid-line matching tests in pmc_interface_mod changed to round-off-error tolerant forms
Property svn:keywords set to `Id`
File size: 220.9 KB

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