Home

Context Navigation

pmc_interface_mod.f90 @ 2240

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

index error in nesting corrected

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

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

Context Navigation

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

Download in other formats: