Home

Context Navigation

pmc_interface_mod.f90 @ 2232

Last change on this file since 2232 was 2232, checked in by suehring, 8 years ago

Adjustments according new topography and surface-modelling concept implemented

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

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

Context Navigation

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

Download in other formats: