Home

Context Navigation

pmc_interface_mod.f90 @ 2223

Last change on this file since 2223 was 2220, checked in by hellstea, 7 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: 175.3 KB

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

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

Context Navigation

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

Download in other formats: