Home

Context Navigation

pmc_interface_mod.f90 @ 2713

Last change on this file since 2713 was 2701, checked in by suehring, 7 years ago

changes from last commit documented

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

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

Context Navigation

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

Download in other formats: