Home

Context Navigation

pmc_interface_mod.f90 @ 2005

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

last commit documented

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/forwind/SOURCE/pmc_interface_mod.f90	1564-1913
/palm/trunk/SOURCE/pmc_interface_mod.f90	merged	eligible
/palm/branches/fricke/SOURCE/pmc_interface_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/pmc_interface_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/pmc_interface_mod.f90	296-409
/palm/branches/suehring/pmc_interface_mod.f90	423-666

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

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

Context Navigation

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

Download in other formats: