Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |