Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |