Home

Context Navigation

pmc_interface_mod.f90 @ 2241

Last change on this file since 2241 was 2241, checked in by hellstea, 8 years ago

last commit documented

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/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: 220.7 KB

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

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

Context Navigation

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

Download in other formats: