Home

Context Navigation

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

Last change on this file since 3792 was 3792, checked in by hellstea, 5 years ago
Interpolations improved. Large number of obsolete subroutines removed. All unused variables removed.
Property svn:keywords set to `Id`
File size: 170.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |