Home

Context Navigation

pmc_interface_mod.f90 @ 2696

Last change on this file since 2696 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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