Home

Context Navigation

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

Last change on this file since 3806 was 3804, checked in by hellstea, 6 years ago
Last commit documented
Property svn:keywords set to `Id`
File size: 170.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |