Home

Context Navigation

pmc_interface_mod.f90 @ 2027

Last change on this file since 2027 was 2020, checked in by hellstea, 8 years ago

last commit documented

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

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

File size: 161.7 KB

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

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

Context Navigation

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

Download in other formats: