Home

Context Navigation

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

Last change on this file since 4413 was 4413, checked in by hellstea, 5 years ago
All the USE-statements within subroutines moved to the module declaration section
Property svn:keywords set to `Id`
File size: 223.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |