Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |