Home

Context Navigation

pmc_interface_mod.f90 @ 2214

Last change on this file since 2214 was 2177, checked in by hellstea, 7 years ago

Minor changes in nesting

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: 174.7 KB

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

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

Context Navigation

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

Download in other formats: