Home

Context Navigation

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

Last change on this file since 4273 was 4273, checked in by monakurppa, 6 years ago
Add logical switched nesting_chem and nesting_offline_chem
Property svn:keywords set to `Id`
File size: 221.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |