Home

Context Navigation

pmc_interface_mod.f90 @ 2903

Last change on this file since 2903 was 2903, checked in by hellstea, 7 years ago

Nesting-related calls to pmci_ensure_nest_mass_conservation and pres after the nest initialization are removed as they may create unwanted initial perturbation in some cases

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/branches/palm4u/SOURCE/pmc_interface_mod.f90	2540-2692
/palm/trunk/SOURCE/pmc_interface_mod.f90	merged	eligible
/palm/branches/fricke/SOURCE/pmc_interface_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90	296-409
/palm/branches/suehring/pmc_interface_mod.f90	423-666

File size: 214.0 KB

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

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

Context Navigation

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

Download in other formats: