Home

Context Navigation

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

Last change on this file since 3803 was 3803, checked in by hellstea, 5 years ago
Nesting anterpolation domain height reduced
Property svn:keywords set to `Id`
File size: 169.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |